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A  Boundary  Element  Technique  in  Transonic  Flow  /I 

Yang  Zuosheng 

(Nanjing  Aeronautical  Institute) 

I.  Introduction 

The  name  "boundary  element"  originates  from  C.A.  Brebbla^^^. 

The  characteristic  of  the  boundary  element  technique  Is  that  all 

the  nodal  unknowns  are  located  on  the  surface  boundary. 

Consequently,  the  numerical  calculation  Is  simplified  and  the 

required  computer  capacity  Is  reduced.  In  this  work,  this 

technique  was  extended  to  the  nonlinear  transonic  flow  about  a 

three-dimensional  wing.  A  weighted  residual  formula  suited  for 

the  transonic  range  was  established  and  directly  applied  to  the 

full  velocity  potential  equation.  In  order  to  ensure  the 

irreversibility  of  shock  vjaves  and  to  facilitate  the  solution  to 

converge  In  the  supercritical  flow  region,  we  added  an  artificial 

viscosity  term  to  the  full  velocity  potential  equation. 

Furthermore,  the  flow  field  was  divided  Into  small  subdomains  and  the 

Green  theorem  was  applied  to  each  element.  If  the  selected 

1 

interpolation  function  has  C  continuity  with  respect  to  the 
velocity  potential  d  for  every  element,  then  the  surface 
Integrals  of  two  neighboring  elements  cancel  each  other  at  the 
Interface  so  that  the  remaining  surface  Integral  (consequently 
the  nodal  unknown)  Is  located  on  the  wing  boundary.  Thus,  the 
boundary  element  Integral  equation  was  derived. 

II..  Full  Velocity  Potential  Equation  with  Artificial  Viscosity 
The  three-dimensional  continuity  equation  is: 

d(pU)  d(pn  d(pfr) 

dx  dy  di  ^  ^  J 

where  the  x-axls  is  the  flow  direction,  y-axis  points  at  the 
right  wing  and  x,  y  and  z  form  an  orthogonal  right-handed 
coordinate.  The  full  velocity  potential  equation  for  a 
nonviscous  potential  flow  can  be  derived  from  equation  (1): 
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-WU^„-WV^.,+  {a^-W’)<f,..  J=  0 

where  P  is  density  and  a  is  speed  of  sound.  U  =  ^  ,  V  =  and  W 

X  y 

s  d  represent  the  three  velocity  components  along  x,  y  and  z 
z 

axes,  respectively. 

In  the  supercritical  flow  region,  the  following  artificial 
compressibility  (viscosity)  was  added  to  prevent  an  expansion 
shock  wave  solution  and  to  eliminate  any  shock  wave  discontinuity 
by  referring  to  equation  (1): 


III "  a  ]  ^  IfH 


where 


M^i  =  AX,  m2  ®  ^y»  »*3 


=  Az.  A  is  a  switch  function. 


A  —  max 


When  the  local  flow  Mach  number  M  is  less  than  ,  A  is  zero.  In 
this  work,  »  0.8.  dp/ax, dp/ay  and  ap/az  represent  the  partial 
differentiation  of  p  in  the  x,  y  and  z  direction,  respectively. 

If  the  artificial  compressibility  term  (3)  is  superimposed  onto 
the  continuity  equation  (1),  the  local  density  will  be  replaced 
by  the  value  upstream  when  the  local  Mach  number  exceeds  M  in 

w 

order  to  permit  receiving  Information  from  upstream.  After 
expanding  equation  (3)  and  adding  it  to  equation  (2),  the 
artificially  viscous  full  velocity  potential  equation  can  be 
expressed  as  (in  tensor  symbols  and  abbreviated  format): 


^  I  ti  —  0  i  l,2f3 


where 


The  repeated  subscripts  indicate  summation  in  the  subscript 
range.  The  subscript  represents  differentiation  0  (when 

i  =  j)  and  1  (when  i  *  j).  ^2  ^3  represents  the  d 

values  at  (x-ax,  y,  z),  (x,  y-Ay,  z)  and  (x,  y,  z-Az), 
respectively. 


III.  Boundary  Element  Integration  Equation 
Equation  (5)  can  be  rewritten  as: 

=  i),  i~6.  ,K,j,  ,=o 

«,/  =  l,2,3  (7) 

Using  a  finite  element  method,  we  divided  the  n  region  of  the 
flow  field  into  M  elements.  If  ^  is  replaced  by  its 
approximation  4  in  each  element,  then  the  weighted  residual 
expression  for  equation  (7)  is: 

i..K„.,GdQ^o  (8) 


For  simplicity,  the  superscript  is  omitted  in  the  following. 

G  is  a  weighted  function  yet  to  be  determined.  We  applied  the  Green 
theorem  to  each  element  0^  in  the  first  volume  integral  on  the 
left  side  of  equation  (8).  If  the  interpolation  function  chosen 

“I 

for  4  has  C  continuity,  i.e.,  4  and  its  first  order  derivative  is 
continuous  across  the  boundary,  then  the  boundary  surface 
Integrals  of  neighboring  elements  will  cancel  out.  After 
applying  Green  theorem  twice  to  the  first  integral  on  the  left 
side  of  equation  (8),  the  surface  Integral  is  located  on  the 
boundary  S.  Thus,  equation  (8)  becomes:  (with  respect  to  wings, 
including  external  surface  and  its  rear  vortex  surface): 


Gf^.,4..nidS-^^  4K,iG,,n.ds'\  +  '^^^^K„G.,.da+ 
+  4G.  iK ,i,  .dQ-'^i  4.  ,K,i,  iGdQ^Q 

a*!'®# 


where  is  the  area  of  the  element  e  on  the  boundary  S.  If  the 

element  0.  is  divided  small  enough,  then  K. ,  can  be  assumed  to  be 
®  th 

a  constant  in  each  element.  In  the  n  iteration,  its  value  is 


3 


equal  to  that  in  the  middle  of  the  element,  ,  in  the 

(n-1)^^  iteration.  Thus,  j=  0  in  each  element.  It 

is  a  step  function.  Its  derivative  K, .  .  or  K. .  .  is  a  pulse 

1  J  j  1  ^  J  ?  J 

function.  In  order  to  satisfy  the  continuity  requirement  for 

across  the  boundary,  we  assumed  that  the  variation  of 
is"^ linear  in  the  vicinity  of  the  boundary.  Thus,  the  volume'' 
integral  in  equation  (9)  is  reduced  to  the  following  surface 
integral  along  the  boundary: 


0G-.  rK.,.  .dQ=^yG. .  AA:,..d5,+  + 

4>.^K.,.,GdQ=^\G4>.  .^K„dS,  +  ^\G4>.  ,\K„dS^  + 

+  udSt 

»•!  ^ 


(10) 


»  =1,2,3 


(11) 


where  P,  Q,  and  R  are  the  local  coordinates  for  the  element.  dSp 

,  dSq  and  dS^  represent  the  differential  areas  on  the  boundary 

when  P  =  const.,  Q  =  const.,  and  R  =  const.,  respectively.  J,  K 

and  L  represent  the  total  number  of  boundaries  on  the  planes 

where  P  =  const.,  Q  =  const.,  and  R  =  const.,  respectively. 

Hence,  equation  (9)  can  be  written  as  the  following  in  the  n^^ 

iteration,  ^rr  f  i  JLr 

2  +  ^K\y'>G.udQ-^ 

••1  ^  J  1  J®, 


•  -1^  fiJ 

-G4>,  ,)'•-*  ’d5<,  +  2jAiC;:-‘  .-G4.  .)'•->  'dS,  =  o 


(12) 


«  ,  i  -1,  2,  3 

If  K. .  is  a  positive  definite  (which  is  always  satisfied  in 
^  J  * 

a  subcrltical  flow) ,  we  obtain  the  fundamental  solution  G  »  to  the 

following  equation  as  the  weighted  function  G: 


is  constant,  the  basic  solution  to  this  equation 


When  K 


vn-  1 ) 

ij 


(13)' 

where  6  is  a  pulse  function. 


i,;  =  l,2,3  (-14) 

is  the  inverse  matrix  of  K . .  and  |  is  the  determinant  of 

ii 

the  matrix  K  .  x^^  and  Xj^'  represent  the  coordinate  vectors  of  a 

fixed  point  and  a  moving  point,  resectively.  When  a  fixed  point 
and  a  moving  point  x^  arc  outside  an  element  then  x^^  ^  ^i  * 

The  third  integral  in  equation  (12)  is  zero.  When  the  fixed 
point  x^  is  inside  an  element  Q^,  based  on  the  characteristic  of 
a  pulse  function,  the  volume  Integral  of  the  element  is  -4lld(Xj^). 

In  the  supercritical  flow  region,  is  no  longer  positive 
definite;  however,  in  equation  (13)  must  be  positive 

definite.  Therefore,  the  third  term  in  equation  (12)  must  be 
rewritten  if  the  neighboring  upstream  element  is  still 
supercritical : 


(ffl=  2 (  <^K'.y"G.udo+H 


(15) 


where 


t  ) 


-R 


(  1  -  1  ) 
i  / 


)G.  i.dQ 


(16) 


Here,  represents  the  value  at  the  first  subcritical 

element  upstream  from  the  supercritical  element.  M  is  the  total 

s 

number  of  supercritical  elements.  Because  M  is  finite,  the 

s 

workload  to  calculate  the  volume  integral  H  is  also  finite. 

Thus,  the  following  boundary  element  integral  for  velocity 
potential  d  can  be  obtained  from  equation  (12): 


where 


=  <f>K\y'  ‘G:,nds]+B-H  (17) 

•  •  I  ^  «  •  I  ^ 


(18) 


AKj^j  represents  the  difference  between  two  neighboring  K^j 
values.  By  respectively  differentiating  equation  (17)  with 
respect  to  x,  y  and  z,  the  corresponding  boundary  conditions  are 


q  .  n  =  0  (on  wing  surface) 

ACp  =  0  (on  tail  vortex  surface) 


(19) 


^  X  (at  infinity) 

where  ACp  is  the  difference  of  pressure  indices  between  the  top 
and  bottom  surface,  and  q  is  the  full  velocity. 

In  Che  transonic  small  perturbation  field  for  a  thin  wing, 
the  boundary  element  integral  equation  (17)  can  be  simplified  as: 


Differentiating  equation  (20)  with  respect  to  x,  we  get 

Differentiating  equation  (20)  with  respect  to  y,  we  get 

where  Xp  represents  the  x  coordinate  on  the  plane  where  P  = 
const.  Differentiating  equation  (20)  with  respect  to  z,  we  get 


(21)  /5 


A:iv(x, 


ill) 


6 


(23) 


dH 

dz 


For  a  small  perturbation  flow  field,  we  obtain  the  following  from 
equation  (6): 


A:.,  =  l-A/;-[3-(2-y)A/i]Afiu(x,y,r)+>4C3-(2-v)^Wi]x 
Af  j/,2)  —  j4[3  —  (2  — v)Af  i]Afiii(jc  — Ajc,y,2) 

K..^  =  Ku  =  l 

K,,  =  K„  =  K,,  =  K,,  =  Ku  =  K,,  =  Q 


(24) 


From  equations  (13)'  and  (14),  we  get 

>  =  l/v/u-i>=  +  A: ’L(y-»7)^  +  (r-i)'] 
In  equations  (20)  -  (23): 


4?)<i,i7,0)  =q>ii,r}r0*'>-<P(^>n,0.) 


(25) 

(26) 


9  is  the  small  perturbation  velocity  potential,  u,  v  and  w  are 
the  small  perturbation  velocity  components  along  x,  y,  and  z, 
respectively,  y  is  the  specific  heat  (y  =  1.4  for  air).  S'^  is 
the  area  of  an  element  e  on  the  z  =  0  plane  where  the  boundary 
surface  S  is  considered  to  be  located  in  approximation.  In  order 
to  limit  the  integration  to  the  wing  region  S^,  we  substituted 
equation  (25)  into  (20).  The  first  integral  on  the  right  side  of 
equation  (20)  was  integrated  by  part: 


“I  )s,,(y-i7)‘  +  ^-  L 


(iS- 


Similarly 


dC 


dSfi-H 


/6 


A«(s,  7, 0) 


(/S- 


•fl. 


ai 

x-i 


dH 

dx 

(28) 


i-vy+2‘2  ] 


y(x-iy‘  +  Klr'‘  C(y-7)'+2^] 

t  ii..  t):::: 


X  ^5p+ 


dy 


(29) 


-  2{(> +'4^)S,,_ X 

’'[‘'^v'(x-iv+A!,'!-~C(li-«‘+»‘3 


X  dS  r  + 


dH 

dz 


(30) 


When  the  fixed  point  (x,  y,  z)  is  located  on  the  wing 
surface  (z  =  0+) ,  the  following  can  be  derived  from  equations 

(28)  and  (30): 


,  ,  dC*” 

4JtCu(x,y,0*) +u(x,y,0-)3=  -  22j\  ^ 


A^dS-rD  +  E 
di 


8 


where 


4«Cw(jf,y,0*)+w<x,y,0.)]=  litn2Sf  [14. 

*-•»  +x'  i 

■^v'u-|)*+A‘.'r~’C(y-7)’+^‘3  (32) 


£■  =  2 


n  c 

dH 
dx 
/ 


G.2^P. 

Oi 


(33)  /7 


For  a  thin  wing,  the  boundary  conditions  (19)  can  be  simplified 
as : 


w  =  dffidz  =  dz'/dx'  -a  (on  wing  surface) 

ACp  =  2  Au  =  0  (on  tail  vortex  surface)  (34) 

V9  =  0  (  ^7z*  *•) 

where  z '  is  the  vertical  coordinate  of  the  lower  wing  surface  in 
the  body  axis  coordinate  x',  y* ,  z*  and  a  is  the  attack  angle. 

The  first  term  on  the  right  hand  side  of  equation  (31)  is  similar 
to  the  linearized  small  perturbation  thickness  integral; 
therefore,  the  result  of  the  linearized  equation  can  be  directly  used. 
Equation  (32)  can  be  considered  as  the  expression  for  a 
linearized  small  perturbation  lift  problem  when  the  effective 
washdown  velocity  at  the  control  point  is  w(x,y,0+)  +  w(x,y,0— )/2 
+  (F+G)/4ll.  Therefore,  once  the  values  of  F  and  G  are 
determined,  it  is  possible  to  solve  equation  (32)  by  familiar 
numerical  methods  for  the  linearized  lift  problem  (such  as 
nuclear  function  method  and  vortex  lattice  method)^^^  to  obtain 
the  corresponding  Au(x,  y,  0)  distribution  along  the  wing  surface 
and  the  F  and  G  values. 


IV.  Solution 

We  divided  it  along  the  wing  surface  into  N  quadrangular 

w 

elements.  Moreover,  planes  (p  =  const.)  perpendicular  to  the 

•>] 


V  V  V 
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wing  surface  were  made  along  the  chord  direction  of  the  wing 

element.  These  planes  extend  to  a  height  approximately  equal  to 

2  times  the  chord  length.  Many  quadrangular  elements  were 

divided  along  the  plane,  as  shown  in  Figure  1.  Figure  1  only 

shows  one  of  the  P  =  const,  planes  and  the  partitioning  of  its 

1 

elements.  In  order  to  ensure  that  v  has  C  continuity  across 
the  boundary  on  a  P  =  const,  plane,  we  expressed  (p  in  an  element 
approximately  as: 


V  LiV,  -r -i- LtV’i-r  LtW, Li<pi L,^Vi 


(35) 


where  (i  =  1,2, - ,  12)  is  the  interpolation  function.  The 

subscripts  1,  2,  3,  4  of  v,  w  represent  the  values  at 
corresponding  nodal  points  of  the  quadrangle,  as  shown  in  Figure 
2.  In  the  element  coordinate  QR  shown  in  the  figure,  L.  can  be 

[31  1 

expressed  as  : 


p-cooi::}l3 


I  ^\J\r' 

m 


f  AR)  f  AQ)  ,  L,^gx(R)gAQ) 
Lx^fAQ)gx{R),  Lx=f,(Q)fAR) 
Ls=g,{Q)fx(R),  Lx’^fAQygAR) 
L,=^f,{R)fAQ),  Lx  =  gx(R)g,{Q) 

L,  =  f,{Q)g,(R),  Ln  =  f,(Q)f,(.R) 
•4u  =  5i(Q)/,(i?),  Z,u  =  /i  (Q)Si(i?) 


1* 

lp» 

(-1,1) 

4 

»  w. 

o 

1 

* 

.  V, 

(-1,-1)  (i,-i> 


(36) 


Figure  1 

1.  p  s  const,  plane 

2.  wing  surface 


Figure  2 


The  functions  are 


The  variable  e  represents  R  or  Q.  The  overall  coordinate  of  a 
point  In  an  element  can  be  expressed  by  the  overall  coordinates 
of  Its  four  nodal  points: 

y  =  (n-Q)y,4. 

J  (37) 

^  ~~(1--R)a-Q)z,  +-J-(i-i?)a+Q)«,+-L<i  +  R)(i+Q),,4. 

The  value  of  u  In  each  element  on  a  P  =  const,  plane  can  be 
approximated  by  the  u  value  at  the  center  of  the  element.  The 
solution  finding  process  Is: 

1.  First,  we  assume  the  u(x,  y,  0^)  values  at  the  center  of 
each  element  on  the  wing  surface,  the  values  of  v  and  w  at  the 
nodal  points  of  each  element  on  F  »  const,  planes,  and  the  values 
of  u  at  the  centers  of  elements.  They  are  used  as  Initial 
values. 

2.  By  substituting  equation  (35)  Into  (33),  the  values  of  D,  /9 
E,  F,  and  G  are  obtained  from  the  given  Initial  values  after 
transforming  the  Integration  variables  to  R  and  Q. 

3.  Using  the  vortex  lattice  method^^^,  the  value  of  u(x,  y, 

0)  at  the  center  of  each  element  on  the  wing  surface  can  be 
obtained  from  equation  (32)  based  on  the  boundary  condition  (34) 
and  calculated  values  of  F  and  G. 


)  1 


4.  The  calculated  values  of  D  and  E,  as  well  as  the  value  of 
u(x,  y,  0-)»  u(x,  y,  0+)-au,  are  substituted  Into  equation  (31). 
Furthermore,  equation  (31)  is  satisfied  at  the  center  of  each 
element  on  the  wing  surface  so  that  u(x,  y,  0  )  can  be  obtained 
at  the  center  of  each  element  on  the  wing  surface. 

5.  According  to  equations  (27),  (29)  and  (30),  the  values  of 
e,v  and  w  are  calculated  at  each  nodal  point  on  P  =  const, 
planes.  The  value  of  u  at  the  center  of  each  element  on  a  P  = 
const,  plane  is  calculated  based  on  equation  (28). 

6.  Steps  2,  3  and  4  are  repeated  until  the  difference  of  two 
consecutive  values  of  u(x,  y,  0+)  at  the  center  of  each  wing 
surface  element  is  less  than  the  specified  value. 

In  order  to  allow  the  solution  to  converge  in  the 
supercritical  region,  the  contribution  from  the  supersonic  point 
(x,  y,  0-)  should  equal  zero. 

V.  Examples 

The  exaiqples  shown  in  Figures  3-5  are  identical  to  those  in 
reference  C4].  The  number  of  elements  on  the  wing  surface  is 
also  equal  to  that  in  reference  C4].  In  the  P  «  const,  plane,  it 
was  partitioned  into  five  regions  along  the  height  according  to 
the  law  of  tangent  (from  0“  -  63.4").  The  partitioning  in  the 
span  direction  is  similar  to  that  on  the  wing  surface.  From 
these  figures,  the  results  in  this  work  agree  well  with  those  in 
references  C4]  and  C5].  However,  when  there  is  an  attack  angle, 
the  front  fringe  result  of  this  work  is  between  those  of 
references  C4]  and  [5].  Figure  6  shows  the  calculated  pressure 
distributions  along  the  :hord  direction  at  three  span  positions 
of  a  dull  leading  edge  rectangular  wing  by  using  this  method  and 
the  difference  method  used  in  reference  [6].  One  can  see  that 
they  are  in  good  agreement. 
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Figure  3.  Rectangular  Wing  with  Aspect  Ratio  4  (double 
arc  wing,  relative  thickness  6%) ,  M  «  0.908 
a  *  0  • 


1.  (wing  tip) 

2.  reference  [5] 

3.  reference  C4] 

4.  this  work 

5.  wing  root 


30  Equi-chord  Length  Sweptback  Wing  with  Aspect 
Ratio  4  (double  arc  wing,  relative  thickness  6X) , 
«  0.908,  a  «  0 

1 .  wing  tip 

2.  reference  [5] 

3.  reference  [4] 

4.  this  work 

5.  wing  root 


Figure  5.  Rectangular  Wing  with  Aspect  Ratio  2  (double  arc 
wing,  relative  thickness  6%) 

1 .  wing  tip 

2.  reference  C5] 

3.  wing  root 

4.  upper  , 

lower  ‘  wing  surface  (ref. [4]) 

5 .  upper , 

lower'  wing  surface  (this  method) 


IS 


Figure  6.  Rectangular  Wing  with  Aspect  Ratio  5  (Model 
NACA64A006 ) 

1 .  this  method 
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A  BOUNDARY  ELEMENT  TECHNIQUE  IN 
TRANSONIC  FLOW 


Yang  Zuosheng 
(Nanjing  Aerenautical  Institute) 

Abstract 

The  boundary  element  technique  is  now  extended  to  study  the  nonlinear 
transonic  flow  about  three  dimensional  wings.  A  weighted  residual  formu¬ 
lation,  valid  throughout  a  Mach  number  range  including  transonic,  is 
developed  and  applied  directly  both  to  the  conviscous  full  velocity  poten¬ 
tial  equation  and  transonic  small  perturbation  equation.  In  order  to 
ensure  the  irreversible  character  of  shock  wave  and  to  make  the  solutions 
stable  and  converged  in  supercritical  region,  an  artificial  viscosity  term 
is  introduced.  We  partition  the  flow  domain  into  a  number  of  small 
elements  and  apply  the  Green  theorem  to  each  element.  The  boundary 
integral  equations  are  obtained  by  using  an  interpolation  function  which 
is  C*  continuous  for  velocity  potential  and  finally  solved  by  means  of 
finite  element  collocation  method. 


The  Exploration  of  the  Spatial  Oscillations  in  Finite 
Difference  Solutions  for  Navier-Stokes  Shocks 

Zhang  Manx in 

(China  Aerodynamic  Research  and  Development  Center) 

Abstract 

In  this  paper,  the  cause  for  the  oscillations  in  the 
upstream  and  downstream  difference  solutions  was  investigated. 

The  study  showed  that  adding  a  second  order  diffusion  term  to  the 
Navier-Stokes  equation  could  smooth  out  the  shock  wave.  However, 
the  third  order  dispersion  term  and  fourth  order  diffusion  term 
could  cause  oscillations  in  the  upstream  and  downstream  solutions 
under  certain  given  conditions.  If  the  second  or  third  order 
difference  method  is  used  to  solve  the  NS  shock  wave  equation, 
the  solution  oscillates  because  of  the  dispersion  and  diffusion 
terms . 

I.  Introduction 

When  a  difference  method  is  used  to  solve  a  shock  wave 

motion,  if  the  difference  method  is  of  the  first  order  of 

accuracy  (or  with  an  added  second  order  diffusion  term),  the 

shock  wave  is  smeared.  If  the  difference  method  is  of  higher 

orders  of  accuracy  such  as  second  or  third  order, oscillations 

[ 1 ] [ 21 

frequently  appear  upstream  and  downstream  .  In  reality, 
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there  is  no  oscillation  .  This  oscillation  still  exists  when 
the  computation  stabilizes.  It  is  significant  to  explain  the 
cause  of  this  oscillation. 

In  this  paper,  an  attempt  was  made  to  perform  an 
enlightening  analysis.  We  believe  that  because  the  step  lengths 
of  spatial  lattice  and  time  are  not  zero  in  the  difference 
equation,  it  is  approximately  correct  to  use  the  difference 
quotient  to  replace  the  derivative.  Therefore,  there  is  a  finite 
difference  between  the  difference  equation  and  differential 
equation.  For  example,  the  model  equation  for  the  initial  value 


problem  is 


du  ,  _  d'u 

dx-'’  dx^  (1.1) 

In  the  corresponding  difference  equation,  when  the  difference  is 
expressed  in  a  Taylor  series  the  equivalent  equation  is: 


du  .  du  d^u  , 


S 


V. 


dx' 


(1.2) 


This  equation  is  called  the  correction  formula  for  (1.1)  where  K 
is  the  accuracy  of  the  difference  lattice  and  v  is  a  coefficient 
related  to  time  and  spatial  steps.  For  compatible  difference 
lattices,  when  the  time  and  spatial  steps  approach  zero,  ■*  0. 
Comparing  (1.1)  with  (1.2),  the  difference  between  differential 
and  difference  equation  is: 


e. 


(1.3) 


This  difference  not  only  affects  the  stability  of  the  difference 
equation  but  also  causes  the  difference  between  the  solution  of 
the  difference  equation  and  that  of  the  differential  equation. 


According  to  this  understanding,  we  used  the  small  perturbation 
method  to  analyze  the  shock  wave  flow  upstream  and  downstream. 

We  studied  the  effect  of  the  second,  third,  and  fourth  order 
terms  on  the  right  hand  side  of  the  correction  equation  to 
illustrate  the  correlation  between  the  difference  format  and  the 
oscillation  in  the  solution.  Some  meaningful  conclusions  were 
also  provided. 

II.  Effect  of  Second,  Third,  and  Fourth  Order  Derivative 
Terms  in  the  NS  Equation  of  Shock  Wave  Motion 
1.  Starting  Equation 

In  order  to  simulate  the  effect  of  ej^  in  the  correction 
formula,  let  us  study  a  one-dimensional  normal  shock  wave  motion 
(see  Figure  1).  The  basic  equation  set  is: 
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/13 


(2.1) 


Here,  p  =  P  ^  the  density  and  the  velocity  of  the 

gas,  t  is  the  time,  x  represents  th  coordinate  (see  Figure  1),  y 
is  the  adiabatic  index,  v  =  (4/3)p.  u  is  the  viscosity  of  the 
gas  which  is  calculated  according  to: 


where 


(A) 


The  subscript  "*  ”  represents  the  value  at  x 
number  and  n  is  the  exponent.  In  equation  (2.1) 


(B) 

M  is  the  Mach 


(C) 


is  the  coefficient  of  the  added  derivative  terms  which  is 

assumed  to  be  a  constant.  It  should  be  pointed  out  that  when  v^ 

's  (i  =  2,3 - )  are  zero,  equation  (2.1)  is  an  accurate  NS 

equation  accurately  describing  the  positive  shock  wave  motion  in 

[3] 

normal  conditions  .  Therefore,  the  steady  state  solution  to 
equation  (2,1)  is  the  correct  shock  wave.  When  v^^'s  are  not  zero 
and  assume  certain  values  related  to  the  lattice  step,  equation 
(2.1)  simulated  the  correction  formula  with  a  specific  difference 
lattice.  Hence,  we  can  study  the  effect  of  second  and  third 
order  terms  by  solving  equation  (2.1). 
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Figure  1.  Velocity  Distribution  of  One-dimensional  Positive 
Shock  Wave 


1 .  real  viscous  shock  wave 

2.  higher  order  difference  calculation 


2.  Effect  of  Second,  Third  and  Fourth  Order  Terms 
If  the  distribution  of  any  physical  quantity  of  a  nonviscous 
positive  shock  wave  is  used  as  the  initial  value,  it  has  already 
been  proven  that  (2.1)  has  a  steady  state  smooth  shock  wave 


solution  when  v„  =  v^  =  v,  0^33.  its  shock  wave  zone  is 

,  2  3  4 

very  thin  (assuming  the  center  is  x  =  0).  Outside  the  shock  wave 


region,  the  physical  quantities  rapidly  approach  the  upstream  and 


downstream  values  of  a  nonviscous  shock  wave.  If  v  ,  v  , -  are 

2  3 

not  zero,  although  (2.1)  still  approaches  its  nonviscous  value 


far  away  from  the  shock  wave  region  yet  some  oscillation  emerges 
near  the  shock  wave  (See  Figure  1).  In  the  following,  the 

V.  and  the  oscillation  is  analyzed. 


correlation  between  v  ,  v  , 


In  reality,  if  we  only  consider  the  presence  of  v^,  v^  and  v 


then  (2.1)  gives  the  following  when  the  solution  is  stabilized: 


(v 


Obviously  u  =  1  and  u  =  U2  are  the  upstream  and  downstream  values 
of  the  nonviscous  shock  wave,  respectively. 

Let  us  study  the  nature  of  the  solution  to  (2.2)  near  u  =  1 


and  u  =  U' 


Because  u  =  1  when  x-  -®  and  u  =  U2  when  x-*  ,  this 


investigation  can  provide  the  characteristics  of  the  solution 
upstream  and  downstream. 

Let  u  =  1  +  u'  (upstream)  (2.5) 

u  =  U2  +  u'  (downstream) 

Assuming  |u'|<<|  in  the  upstream  and  |u‘|<<U2  downstream, 

equation  (2.5)  is  substituted  into  (2.2).  After  omitting  higher 

2 

order  small  terms  beyond  u'  ,  we  get 


du' 


d*tt'  d‘u'  f  , 


(upstream) 
(downstream)  (2.6) 


where 


-U2)>0 


(2.7) 


v^^  and  V2^  are  the  v  values  at  u  =  1  and  U2 ,  respectively. 
Obviously,  (2.6)  is  a  linear  equation  whose  solution  can  easily 
be  determined  by  the  following  characteristic  equations: 


(upstream) 
(downstream)  (2.8) 

The  discussion  is  carried  out  in  the  following  cases: 

(1)  w-j  and  ^2  greater  than  zero  and  v^  =  =  0.  The 

characteristic  root  of  (2.8)  is 


A  =  — 

/It 

X - ^ 

/'l 


(upstream) 


(downstream) 


The  general  solution  to  (2.6)  is: 


Or,  one  can  get  from  (2.5)  that 


ii'  =  • 


(upstream) 


(downstream) 


a  *1  +Aei‘,  * 


(upstream) 

a  u,+Be  (downstream)  (2.9) 
(2.9)  shows  that  the  upstream  solution  smoothly  approaches  1 
exponentially  when  x-  -•  .  When  x  -•  ,  the  downstream  solution 
snoothly  approaches  U2  wihtout  any  oscillation. 

(2)  U‘j>0,  “  0 

From  (2.8)  in  the  upstream  if 


(2.10) 


then  both  characteristic  roots  are  positive  real  numbers.  The 
general  solution  to  (2.6)  is: 

H^l  +  A .  «p[(-^  +  - 4v,fc,  )  *  ]  + 

r,  ,  ,  ,  (2.11) 

)*  1 

Here,  and  A2  are  integration  constants  and  the  solution  does 
not  oscillate.  If 


both  characteristic  roots  of  (2.8)  are  complex  numbers. 
Therefore,  the  general  solution  of  (2.6)  is: 


(2.12) 


This  equation  shows  that  the  solution  exhibits  oscillations  in  the 
upstream.  Moreover,  the  amplitude  gradually  increases  with 
increasing  x.  Figure  2a  shows  this  change. 

In  the  downstream  region,  because  all  characteristic  roots 
of  (2.8)  are  real  and  u  -*  U2  when  x  -*•  ,  the  solution  to  (2.6) 
is : 

1 

5  «  S ,  +  ^  ex p{[-^  — ^  (/i!  Ak,v,)  ]  X  I  (2.14) 

One  can  see  that  the  solution  does  not  oscillate. 

If  we  assume  that  ^^<0  and  =  0.  The  solution 

does  not  show  oscillation  upstream  as  one  can  see  by  the  same 
analysis.  In  the  downstrean.  however,  if 


the  solution  shows  osclllatlons(see  Figure  2b). 


O 


Figure  2a  Uo  >0,  v-  >0  and  v,  =  0 


Figure  2b  W2  '^4®  ® 

(3)  v^  >0,  v^  =  =  U2  ®  ® 

In  the  upstream  region,  the  characteristic  roots,  given  by 
(2.8)  are: 


where 


a,=a;‘>  +A“>f 
'A,  =  Ar' 


(2.15) 


(2.16) 


Considering  u  -  1  when  x  -  the  solution  (2.6)  is: 

u-l+A,e*^  ‘cos(A,“*jr)  'sin(A,''*x)  (2.17) 

Here,  and  A2  are  integration  constants.  This  solution 
obviously  oscillates.  Furthermore,  the  amplitude  increases  with 
increasing  x  (see  Figure  3). 

For  the  downstream  region,  the  solution  of  (2.6)  can  be 
obtained  similarly: 


u  =  u  i-r  'cos(A‘  '  sinCA/’’ jc)  (2.18) 


where 


(2.19) 


It  shows  that  the  solution  also  oscillates  downstream.  Moreover, 
the  amplitude  Increases  with  decreasing  x  (See  Figure  3). 

(4)  v^,  V2,  U”!  and  U2  simultaneously  are  present.  /I 

Assuming  v^>  0,  based  on  the  determination  formula  of  the 
third  order  equation  one  can  easily  find  out  that  If 

+  (2.20) 

then  the  solution  of  (2.6)  oscillates  upstream.  In  the 
downstream  area,  if 

(2.21) 

then  the  solution  of  (2.6)  also  oscillates.  From  (2.20)  and 
(2.21)  one  can  see  that  when  V2>  0,  it  0  and  u-j  =  112  ® 
always  oscillates  upstream.  In  the  downstream,  however,  only 
when 

ktvl  4 
vj  ^^27 


it  exhibits  oscillation.  If  v^  <0,  v^  #  0  and  u-j  =  U2  *  0>  there 
arealways  oscillations  downstream.  Upstream,  however,  only  when 


it  exhibits  oscillations (see  Figure  4). 
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III.  Discussion  on  Upstream  and  Downstream  Shock  Wave 
Oscillation  In  Difference  Calculation 
Based  on  the  above  analysis,  one  can  see  that: 

1,  In  the  viscous  shock  wave  propagation  equation,  a  smooth 
solution  can  be  obtained  by  adding  ^2  alone.  The  second  order 
diffusion  term  can  smooth  shock  waves.  If  second  and  third  order 


terms  are  added ,  when  v^  >  0  and 


the  upstream  shock  wave  solution  shows  oscillations  and  the 
downstream  solution  is  smooth.  If  v^CO  and  _  /18 


the  downstream  solution  exhibits  oscillations  and  the  upstream 
solution  Is  smooth.  If  the  fourth  order  term  Is  added.  Its 
effect  Is  to  simultaneously  cause  or  aggrevate  oscillations  In 
the  upstream  and  downstream. 

2.  Because  the  second  order  lattice  has  a  third  order 
dispersion  term  and  a  fourth  order  diffusion  term  and  V2  -  0,  the 
upstream  and  downstream  shock  wave  exhibits  oscillations  when  v  = 
0  or  very  small.  If  the  effect  of  the  fourth  order  diffusion 
term  is  far  less  than  that  of  the  third  order  dispersion  term, 
then  one  side  of  the  shock  wave  exhibits  larger  oscillations 
depending  on  the  sign  of  the  dispersion  term.  In  order  to 
minimize  oscillations,  we  should  try  to  reduce  the  third  and 
fourth  order  terms.  For  the  explicit  MacCormack  form,  if  the 
Courant  number  approaches  1  and  v^  and  v^  are  small,  the 
oscillation  should  be  small.  It  has  already  been  proven  In  the 
numerical  experiment  In  reference  [1]. 

3.  The  fourth  order  diffusion  term  appears  because  V2  =  v^ 

=  0  in  the  third  order  lattice.  Therefore,  it  is  unavoidable  to 
have  oscillations  In  the  upstream  and  downstream  shock  wave  when 
V  *  0  or  is  very  small.  However,  it  is  generally  smaller  than 
the  second  order  oscillation.  To  further  minimize  osclllatlonsi 
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we  should  reduce  one  fourth  order  diffusion  term.  It  is  expected 
that  the  fourth  order  lattice  only  has  negligible  oscillations  in 
the  upstream  and  downstream  because  V2  *  =  0.  The  effect 

of  higher  order  dispersion  and  diffusion  terms  is  very  small. 

4.  In  order  to  verify  the  accuracy  of  these  conclusions,  we 
gathered  some  numerical  shock  wave  experiments  in  the  literature 
based  on  the  following  Burgers  equation: 


Figure  5  is  the  calculated  result  given  in  reference  [4]  using 
the  first  order  difference  lattice.  When  a  =  0.038925,  v  =  o  ind 
X  =  u  =  2.593467.  When  x  =  •,  u  =  1.393784.  One  can  see 
that  the  shock  wave  curve  is  smooth.  Figure  6  shows  the 
calculated  results  given  in  reference  [5]  using  the  explicit 
MacCormack  second  order  lattice.  When  a  *  0,  v  =  10”^  and  x  *  -« 

,  u  *  -1.  When  x  =  •,  u  =  1.  One  can  see  that  oscillation 
occurs  both  upstream  and  downstream.  Figure  7  shows  the  results 
given  in  reference  [4]  using  the  third  order  lattice.  The 
conditions  are  identical  to  those  in  Figure  5.  One  finds  that 
there  are  oscillations  on  either  end  of  the  shock  wave.  However, 
it  is  much  smaller  than  that  using  the  second  order  lattice. 

Figure  8  shows  the  results  given  in  reference  [6]  using  the 
fourth  order  lattice.  When  a  =  -(1/2)  and  x  =  -•,  u  *  1.  When  x 
=  •,  u  =  0.  The  lattice  Reynolds'  number  is  2.5.  One  can  see 
that  there  is  almost  no  oscillation  in  the  shock  wave. 
Correspondingly,  if  the  second  order  lattice  is  used,  oscillations 
will  appear  (see  C6]). 
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Figure  5. 


Burgers  Equation  Shock  Wave  Solution  Given  by 
First  Order  Difference  Lattice[4] 


1.0 

0 

-1.0 

Figure  6.  Burgers  Equation  Shock  Wave  Solution  Given  by 
Second  Order  MacCormack  LatticeCS] 

Similar  conclusions  are  derived  by  numerical  simulation  of  the 

[1  21 

Euler  equation  solution  of  the  shock  wave  ’  . 
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Figure  7.  Burgers  Equation  Shock  Wave  Solution  Given  by 
Third  Order  Lattice[4] 


Figure  8.  Burgers  Equation  Shock  Wave  Solution  Given  by 
Fourth  Order  Lattice[6] 

Results  of  numerical  simulation  completely  confirmed  the 
accuracy  of  the  conclusions.  Therefore,  the  cause  of  oscillations 
in  upstream  and  downstream  shock  wave  solutions  using  a 
difference  method  is  explained. 

During  the  course  of  writing  this  paper,  the  author 
discussed  with  Comrade  Gao  Shuchun  many  times.  He  wishes  to 
express  gratitude  for  his  assistance. 
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THE  EXPLORATION  OF  THE  SPATIAL  OSCILLATIONS 
IN  FINITE  DIFFERENCE  SOLUTIONS  FOR 
NAVIER-STOKES  SHOCKS 


Zhang  Hanxin 

(China  Aerodynamic  Research  and  Development  Centre) 

In  this  paper,  the  spatial  oscillations  in  finite  difference  solutions 
for  Navier-Stokes  shocks  are  explored.  It  is  shown  that  the  second  order 
diffusion  term  added  to  NS  equations  could  smear  the  shock  wave,  damp 
the  oscillations  in  the  vicinity  of  the  shock.  However,  the  third  order 

dispersion  term  and  fourth  order  diffusion  term  added  to  NS  equations  could 

cause  the  oscillations  in  upstream  and  downstream  region  of  the  shock 
under  the  conditions  given  by  this  paper.  Therefore  the  oscillations  in  the 
difference  solutions  with  second  or  third  order  accuracy  could  arise  from 
the  numerical  dispersion  and  diffusion  terms. 
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Calculation  of  Boundary  Layer  Growth  Behind  An  Unsteady 
Expansion  Wave  In  a  Tube 
Wang  Songgao 

(Institute  of  Mechanics,  Chinese  Academy  of  Sciences) 

I.  Introduction 

In  order  to  improve  the  operating  quality  of  wind  tunnels 
and  shock  wave  tubes,  the  study  of  boundary  layer  growth  after  an 
expansion  wave  enters  a  cylindrical  tube  is  of  significance. 
Especially  in  transonic  wind  tunnels,  the  boundary  layer  growth 
in  the  gas  storage  tube  with  time  will  directly  affect  the 
operating  time  and  flow  quality  of  the  tunnel.  To  solve  this 
problem  is  the  basis  of  wind  tunnel  design. 

r  1 1 

E.  Becker '■  ■'  first  calculated  the  boundary  layer  growth  in 
the  gas  storage  tube.  He  used  the  two  element  incompressibility 
assumption  and  the  1/7th  power  velocity  distribution  to  simplify 
the  boundary  layer  momentum  integral.  Then,  the  solution  was 
found  by  using  the  Blasius  skin  friction  law.  Furthermore,  an 
effective  central  expansion  wave  velocity  was  obtained.  Although 
the  effect  of  compressibility  and  heat  transfer  on  p/p^  was 
considered  in  the  boundary  layer  and  the  skin  friction 
coefficient,  the  result  is  still  incompressible.  H.  Ludwieg  C^l 
modified  the  Becker  method  for  the  velocity  cross-section  and 
skin  friction  coefficient  to  obtain  better  results. 

J.C.  Slvells'  work^^^  pointed  out  that  corrections  must  be 
made  in  axisymmetry,  skin  friction  coefficient,  velocity  cross- 
section  and  effective  expansion  wave  origin.  Becker's  effective 
expansion  wave  propagation  speed  could  still  be  used.  In  his 
calculation,  two  modifications  were  made:  one  is  to  modify  the 
result  of  a  flat  plate  and  the  other  is  to  correct  the  origin  of 
the  effective  expansion  wave  according  to  the  start-up  time  of 
the  experiment.  A  numerical  method  was  used  in  the  calculation. 

We  began  directly  from  the  unsteady  axisymmetric  boundary 
layer  momentum  integral  to  find  the  boundary  layer  growth  in  and 
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behind  the  expansion  wave  assuming  that  the  steady  state  boundary 
layer  velocity  crss-section  and  skin  friction  could  be  applied  to 
an  unsteady  boundary  layer,  and  the  axial  pressure  gradient  of 
the  nuclear  flow  in  the  tube  is  negligible.  Finally,  the  formulas 
to  calculate  the  effective  expansion  wave  propagation  speed  and 
boundary  layer  thickness  were  found.  These  formulas  are  in 
algebraic  expressions.  In  the  binary  incompressible  case, 

Becker's  results  are  obtained.  The  computation  process  is 
simple.  It  is  in  good  agreement  with  the  experimental  result. 


For  engineering  design,  modification  of  the  origin  of  the 
effective  expansion  wave  can  be  skipped. 


II.  Basic  Equation  and  Its  Solution 
The  unsteady  boundary  layer  momentum  integral  for  a  tube  is 


where 


dx  u.  dx  ^  dx  ^  r. 


(1) 


(2) 


Here,  p,  u  and  rg  are  the  density,  flow  velocity  and  inner  radius 
of  the  tube.  The  subscript  e  represents  the  free  flow  direction 
of  the  boundary  layer,  c^  is  the  skin  friction  coefficient,  x 
is  the  axial  coordinate  and  its  positive  direction  coincides  with 
the  propagation  direction  of  the  expansion  wave,  y  is  the  radial 
coordinate  which  points  from  the  wall  toward  the  center,  t  is 


From 


M,  dx  uj  dt  M,  \  ox  at  / 


=.-A\  ^=0 

p.u:  Ox 

Equation  (1)  can  be  converted  into: 

dx  ‘  u,  dx  ‘  p;  dx  ■'■p.u.  dt 

Substituting  the  power  law  velocity  distribution  into  equation 
(2),  we  get 
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When  n  -  ®,  the  above  formula  approaches  1.  When  n  =  7~11,  the 
change  caused  by  (T^/T^)  =  0.8 —  1.0  and  (6/rQ)  =0-1  is  around 
47..  The  entire  change  is  approximately  81.  Because  we  assumed 
that  A  is  a  constant,  equation  (3)  can  be  rewritten  as: 


Let 


(4) 

(5) 


The  subscript  4  represents  the  gas  storage  state.  Equation  (4) 
is  converted  into 


Pt.  +  A-  +  A  ^  A.  P.“« 

dx  u.  dx  u.  \  dt  u,  dt  /  2  p.u, 


Equation  (6)  is  the  basic  equation  of  this  discussion. 

Generally,  an  unsteady  expansion  wave  is  a  wave  series  of 
finite  thickness.  The  boundary  layer  growth  in  the  wave  and 
behind  the  wave  should  be  separately  discussed.  The  two 
solutions  are  linked  together  at  the  tail  of  the  wave.  Figure  1 
illustrates  this  situation.  In  the  following,  we  will  discuss 


in  two  separate  situations  and  then  combine  them  together. 


Figure  1.  Schematic  Diagram  of  Boundary  Layer  Growth 

1.  composite  solution 

2.  connecting  point 

(1)  Central  Wave 

When  the  thickness  of  an  unsteady  expansion  wave  is  zero, 
the  nuclear  flow  parameters  p^  and  u^  in  the  tube  behind  the  wave 
are  invariant.  Equation  (6)  can  be  simplified  as: 


dtf  1 

dx  U.  dt  2  (7) 


According  to  the  theory  of  partial  differential  equation,  we  can 
obtain  the  characteristic  equations  for  equation  (7): 


dx 

ds 


1 


ds 


(8) 


Here,  V  represents  the  propagating  speed  of  a  central  expansion 
wave.  From  the  above  formula,  we  can  obtain  an  equivalent 
differential  equation  for  equation  (7): 


X  =  x-yt 


After  ascertaining  =  c^Ce),  equation  (9)  can  be  solved. 
(2)  Inside  the  Wave 


Let 


(10) 


Here,  a  is  the  speed  of  sound.  The  subscript  4  represents  the 
gas  storage  state.  We  get 

^  _  d  _  1  — d 
dx  di  i  dri 
1  4  -d-?)*  d 
a,  dt  i  dr) 

When  y  =  1.4,  the  correlation  between  parameters  inside  an 
unsteady  expansion  wave  and  those  in  front  of  the  wave  is^^^: 


I  (11) 


We  have 


The  viscosity  was  approximated  by 


the  Sutherland  formula 


[ID: 
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Using  the  above  formulas,  equation  (6)  can  be  written  as; 


(15) 


Because  Re^  =  PgU^e  =  4>  p^a^  =  f(^,n),  and  (^,n)  in 

“4  “e 

equation  (15),  the  characteristic  equation  set  for  equation  (15) 
is 
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From  the  two  equations  above,  we  get 
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The  solution  is 


(17) 


The  constant  c  can  be  determined  later. 

From  d^/ds  =  d^/dy  dn/ds  and  equation  (17),  the  equivalent 
normal  differential  equation  for  equation  (15)  is: 


j*  * 


(‘ 


(18) 


After  C£  is  determined,  d-d(n),  i.e.,  e  se(n),  can  be  determined. 

(3)  Composite  Solution 

Strictly  speaking,  the  boundary  layer  growth  process  should 
be  solved  inside  and  outside  the  wave  as  shown  in  Figure  1.  The 
complete  solution  is  obtained  by  combining  the  two  solutions  at 
the  end  of  the  wave.  Obviously,  solving  equation  (18)  is  most  of 


the  work.  We  can  assume  a  composite  solution:  to  find  an 
equivalent  expansion  wave  of  zero  thickness  which  results  in  the 
same  boundary  layer  thickness  at  the  tail  end  of  the  wave  at  the 
same  time  as  that  in  equation  (18).  We  used  such  a  zero 
thickness  expansion  wave  to  find  a  universal  solution.  In  fact, 
it  is  a  problem  to  find  the  propagating  speed  V  of  an  equivalent 
central  wave. 

Assuming  that  the  friction  coefficient  law  is  normal:  /24 


c, =&(/?«,)--  (19) 

k  is  a  constant  and  the  subscript  0  represents  taking  Re  with 
respect  to  e. 

Substituting  it  into  equation  (9),  we  get  the  solution 
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into  it, 

we  get: 


V 
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(21) 


At  the  tail  n  =  n-j  =  Z  j  ,  generally,  M2<  0.30  which 

1+  1^  M3 

means  n-i  <0.34.  n-j  is  a  small  quantity.  After  keeping  the 
lowest  order  term  on  the  right  side  of  equation  (21),  we  get: 


l+m  \  72  * 

The  subscript  1  represents  the  wave  tail. 
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Corresponding  to  equation  (18),  we  get 
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The  solution  is 
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At  the  tail  of  the  wave: 
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From  5/1-n  =  a^t,  we  can  specify  that  c  =  a^t  C  1-(  1-5/6  a) 
Substituting  it  into  the  above  equation  and  keeping  the  lowest 
order  term  in  n-],  we  get: 
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Comparing  (22) to  (25)  we  get: 


i.e. , 
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When  the  flow  Is  binary  and  incompressible  and  n=7,  m=1/4,  a 
=  n+2/n  =  9/7. 


(28) 


This  is  Becker's  result^ Strictly  speaking,  the  propagating 
speed  V  of  the  axisymmetric  effective  central  expansion  wave  is 
also  related  to  the  velocity  cross-section  and  the  law  of  skin 
friction  in  first  order  approximation. 

III.  Results  and  Discussion 

For  convenience,  we  take  the  absolute  value  of  u^.  Then, 

e 

equation  (9)  can  be  rewritten  as; 
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where 
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As  long  as  c^  is  determined,  the  solution  can  be  found.  / 

We  noticed  that  the  value  of  1/2-m  varies  by  about  3X  when  m 

=  1/4  ~ 1/7.  We  also  assumed  that  6  -6  /e  is  a  constant 

P  ^ 

(although  the  not  too  large  variation  of  a  -Ap/B  is  included  in 
the  calculation),  equation  (26)  or  (27)  may  be  considered 
invariant  in  approximation.  Therefore,  the  assumption  that 
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(33) 


is  not  related  to  the  specific  laws  of  velocity  cross-section  and 
skin  friction  in  reference  [3]  is  still  valid  (approximately)  in 
axlsymmetric  conditions.  In  the  normal  range,  the  difference  in 
the  values  of  (31)  and  (33)  is  less  than  107.. 

Equations  (29)  to  (32)  are  our  starting  equations.  For 
simplicity,  equations  (31)  and  (32)  can  be  approximated  by  the 
corresponding  Becker's  result. 

Let  us  take  the  Karman  modified  friction  coefficient 
expression 


_ (0.2*2)  • _ ^xV-1 

(Ig /?e, +1.1696)  +0.3010)  2  ^*/ 


(34) 


and  substitute  it  into  equation  (29)  to  obtain  the  boundary  layer 
thickness  by  integration: 
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Subsequently,  the  expressions  for  9/6,  6  /6  and  6/6  are  found. 
After  reorganization,  the  above  formula  can  be  written  as: 


(36) 
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.If  Tg/X^  =  1,  the  above  formula  can  be  further  simplified. 
Obviously,  the  computation  is  simple. 
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Figures  2  and  3  are  the  calculated  results. 

Figure  2  is  a  comparison  of  the  calculated  results  with 
experimental  values  in  wind  tunnels.  They  agree  very  well. 
Results  calculated  by  Ludwieg  are  also  shown  in  the  figure. 


Figure  2.  Comparison  and  Calculated  Results 

1 .  this  work 

2.  t  (microsecond) 


Figure  3.  Comparison  and  Calculated  Results 

1 .  this  work 

2.  t  (microsecond) 


43 


Figure  3  is  the  comparison  of  our  calculated  results  wth  the 
experimental  wind  tunnel  values  at  Arnold  Engineering  Development 
Center  (AEDC).  The  calculated  results  deviate  from  the 
experimental  values  by  a  certain  moment »  In  the  figure,  the 
modified  result  in  reference  [3]  is  also  given.  Sivells  used  the 
start-up  time  of  the  experiment  to  modify  the  effective  origin. 
This  is  to  add  another  time  correction.  When  the  starting  device 
is  located  downstream  from  the  wind  tunnel,  this  correction  could 
improve  the  agreement  with  experimental  values.  However,  this 
correction  could  not  be  obtained  ahead  of  time  in  engineering. 
Numerically,  this  correction  value  would  not  bring  about  a  large 
error  for  design  purposes.  It  could  be  neglected.  Of  course,  it 
is  ideal  to  determine  this  correction  in  theory. 

In  summary,  it  is  possible  to  theoretically  calculate  the 
boundary  layer  growth  behind  an  -unsteady  expansion  wave  in  a 
cylindrical  tube  by  using  the  composite  solution  of  an  effective 
central  wave  propagating  at  a  speed  V  directly  from  the  unsteady 
axisymmetric  boundary  layer  momentum  integral  with  some 
assumptions.  The  effective  velocity  V  has  been  found.  Within  a 
certain  error  range,  it  is  not  related  to  the  velocity  cross- 
section  and  the  law  of  skin  friction.  Becker's  result  is  a 
special  case  obtained  in  a  binary  incompressible  system  when  n  = 
(1/7)  and  m  =  (1/4).  The  composite  solution  is  an  algebraic 
expression  which  is  easy  to  calculate.  The  results  are  in  good 
agreement  with  the  experiments.  For  engineering  design  purposes, 
the  accuracy  is  sufficient. 
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CALCULATION  OF  THE  BOUNDARY  LAYER  GROWTH 
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Abstract 

The  problem  of  the  boundary  layer  growth  in  the  charge  tube  is  dis¬ 
cussed.  Based  on  E.  Becker's  work  and  J.  C.  Sivells'  modification,  in  the 
axisymmetric  case,  a  theoretical  treatment  of  the  boundary  layer  growth 
in  the  expansion  wave  is  giveni  the  movement  velocity  of  the  equivalant 
unsteady  expansion  wave  of  zero  width  is  derived;  an  analytical  solution 
-s  obtained  and  has  been  reduced  to  an  algebraical  expression.  |  The  result 
contains  various  factors  which  affect  boundary  layer  growth;  axisymmetry, 
velocity  profile,  skin-friction  coefficient  law  and  expansion  wave  thickness 
effect.  The  calculation  is  simple  and  the  results  coincide .  with  experi¬ 
ments. 
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Abstract 

In  this  work,  a  computation  method  for  inviscid  supersonic 
flow  around  a  bent  cone  was  developed  using  a  difference  method. 
In  order  to  overcome  the  difficulty  brought  about  by  the  bent 
body  axis  in  the  calculation,  the  equation  solving  place  is 
transformed  into  an  inclined  plane  in  the  advancing  direction 
through  the  Euler  transformation  of  the  independent  variables. 

In  other  words,  an  "inclined  shift"  or  "parallel  shift"  of  the 
Euler  equation  was  used  to  solve  the  problem.  Numerical  examples 
showed  that  good  results  could  be  obtained  using  this  method. 

I.  Introduction 

This  numerical  method  for  the  characteristics  of  an  inviscid 
supersonic  flow  around  a  bent  cone  was  based  on  a  difference 
method  and  the  shock  wave  capturing  technique  with  the  following 
essential  points.  (1)  In  order  to  overcome  the  difficulty  due  to 
the  bent  cone  axis ,  the  independent  variables  of  the  Euler 
equation  are  transformed  onto  an  inclined  solving  plane  which  is 
the  same  direction  as  the  thrust.  Moreover,  the  angle  of 
inclination  gradually  increases  or  decreases  as  the  confutation 
progresses.  In  other  words,  the  numerical  solution  of  the  Euler 
equation  is  found  by  the  "inclined  advance"  with  a  varying  angle 
of  inclination  or  "parallel  advance"  using  right  inclined  angle. 
Its  advantage  is  that  the  calculation  can  advance  in  a  rational 
direction  according  to  the  specific  shape  of  the  object  to 
simplify  the  calculation  and  to  ensure  the  accuracy.  (2)  The 
second  order  two-step  MacCormack  scheme  is  used.  (3)  The 
computation  area  is  divided  into  two  parts  by  a  transformation 
plane  (see  Figure  1).  The  part  in  front  of  the  plane  uses  the 


front  cone  cylindrical  coordinate  to  solve  the  Euler  equation  by 
"inclined  advance".  The  part  behind  the  plane  uses  the  rear  cone 
cylindrical  coordinate  to  solve  the  Euler  equation  by  "parallel 
advance".  (4)  The  computation  formulas  for  the  external  shock 
wave  and  the  boundary  are  derived  from  the  shock  wave  correlation 
and  the  boundary  compatibility  relations,  respectively.  (5)  The 
flow  characteristics  of  singular  points  on  the  intersect  of  the 
front  and  rear  cone  surface  are  solved  individually.  The 
position  of  internal  shock  wave  is  derived  from  the  pressure 
charge  in  the  flow  field.  No  filtering  and  smoothing  process  is 
included  in  the  calculation.  (6)  A  non-uniform  radial  lattice  is 
used  so  that  the  number  of  meshes  is  Increased  near  the  surface 
without  increasing  the  total  number  of  meshes. 

II.  Basic  Equations  and  Boundary  Conditions 

The  coordinate  system,  attack  angle  a.  and  sideslip  angle /9 

% 

used  in  the  computation  are  defined  in  Figure  1.  Different 
equations  were  used  in  the  "inclined  advance"  and  "parallel 
advance"  areas  for  solving  the  Euler  equation.  In  the  latter 
case,  the  conserved  Euler  equation  is  used.  The  specific 
computation  method  is  shown  in  reference  [2].  In  the  following, 
the  equation  in  the  "inclined  advance"  region  and  its  boundary 
conditions  are  given. 

1.  Points  in  the  Flow  Field 

In  a  cylindrical  coordinate  (z,  r,  ^),  the  aerodynamic 
equation  for  a  steady  invlscid  and  thermally  non-conductive  flow 
is : 
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The  pressure  p,  density  p,  speed  of  sound  a,  and  velocity 
components  u,  v,  w  are  dimensionless.  The  dimension  factors  are 
Pc^,V^  »P»  order  to  change  the  solving  plane  of 

equation  (1)  from  z  =  const  to  an  inclined  plane  ri(z»  r,  = 
const  (with  respect  to  z-axis),  new  variables  n,  C,  and  9  are 
introduced  in  the  following  equation. 

z  =  £{r,ij,9)  ^rj  +  riA  +  Bri^costp 

where 

A  =ctg9t — — (ctgfl*-ctgflj) 

^  Vq 

B  =  -5-^ — (ctgff*-ctg  9,) 

9  7o 

where  nn  and  Qr>  the  initial  values  of  n  and  angle  of  inclination 

^  ^  if  ic 

of  the  front  cone  flow  field  solution,  n  and  0  represent  the  n 
value  and  angle  of  Inclination  of  the  transformation  plane  (see 
Figure  1).  f^(n,9)  and  f2(n,9)  represent  the  external  shock 

wave  boundary  and  object  boundary  [r=f-j(n,9)»  r=f2(n>9)]» 


respectively. 

By  using  the  transformation  correlation  (2),  equation  (1) 
finally  becomes 


^  +»7i  +  (4- 


where 


a/> 

dri  o*« 


w 

au  1 

r 

dt  P 

w 

dv  1 

r 

di>  P 

w  au)  : 

r  df  f. 

I  ^  “  V  a{  +  r  ;J 


^  ’V  -TI,U  +  t)r  «+— JJ,U> 


2  'V  =4.«+  6rV  +  ~i^W 

dij  1 
dz  “E. 


on 

£,=  JE^iJlil  =  U  +  Bn)co^i> 

or 

J-£,=J_.M%^  =-JLu  +  Bn)^.ni> 

r  T  dip  r 

Equation  (3)  is  the  formula  to  calculate  a  point  in  the 
"inclined  advance"  solving  region.  When  the  solving  plane 

“ic 

advances  from  no  to  n  ,  its  angle  of  inclination  with  the  z-axis 
also  gradually  changes  from  0o  to  e*. 

2.  Shock  Wave  and  Characteristic  Compatibility 
According  to  an  analysis  of  the  shock  wave  boundary 
characteristics,  the  only  characteristic  compatibility  in 
calculating  the  shock  wave  point  is  the  Porter  I  family 
compatibility  relation. 


where 


If  y  -  ]  -  <  If  -  ’--If "  7-  '  If  -  ^ 

F  =  ~ paiN  iG  1  +  N ,G i + N ,G ,)  +a^Gi  +(?, 


,V, 

\r»  rW  v'-v\)  sv'i-vl  " 


y  =  y/  M* +  !/*  +  «/' 

V  ,==un„+un„  +  u;n,, 


n„  =— - =====  1 

%/l  +mj  +mj 


where  n  n  and  n  ,  are  the  unit  normal  vector  of  the  shock 
si’  s2 ’  s3 

wave  plane, 


dz  ’  f, 
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In  addition,  the  shock  wave  point  solution  should  also 
satisfy  the  Ranke-Hugoniot  relation: 

y+l  V(y+1)'  +  ’ 

P=-L-.piyi^ 

s—p 

(5) 

where 

S^p,-hp,r:, 

p-j,  p^,  u^,  v^,  and  w^  are  quantities  in  front  of  the  wave. 

In  order  to  facilitate  obtaining  the  numerical  solution,  the 
Ranke-Hugoniot  was  differentiated  with  respect  to  n.  It  was  then 
used  simultaneously  with  the  characteristic  compatibility 
relation  to  derive  the  differential  equation  for  the  shock  wave 
slope  m^ 

dfiii  — 

drj  Lm,  (6) 


where 

L.,  =po  ]^{eK,  +(eK^  +_L  +  - 

-a'f{p.K,-\-r),K,  >7,  K^-eK^  j 

I,,  *po  [(e/,  +-L  +^e/,  +  _L  ri,I ^  N,^{eK, 
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All  the  derivatives  of  flow  parameters  with  respect  to  m-j  and  m2 
can  be  obtained  from  the  Ranke -Hugoniot  relation. 

Equations  (5)  and  (6)  are  the  equations  to  solve  the  shock 
wave  point. 

3.  Surface  Boundary  Conditions  and  Characteristic 
Compatibility  Relations 

From  an  analysis  of  the  surface  boundary  characteristics  we 
know  that  there  are  four  compatibility  relations  in  calculating  a 
surface  point:  three  flow  characteristic  compatibility  relations 
and  the  Porter  II  series  compatibility  relation.  They  are 
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E i^GiU  +  G jV +  GiXj!  ,  E2^Giti +Gifi +  G4fi 
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n^,  n2,  and  n3  are  the  unit  normal  vectors  of  the  surface 
In  addition,  a  surface  point  should  also  satisfy  the 
following  boundary  condition; 


III.  Numerical  Computation 
The  computation  of  a  supersonic  inviscid  flow  is 
mathematically  a  problem  to  solve  the  initial  boundary  value  of  a 
five  variable  quasi>llnear  hypobollc  equation  series.  After  the 
solution  on  the  given  initial  plane  is  obtained,  the  computation 
can  advance.  Specifically  for  a  bent  cone,  this  method  involves: 

(1)  using  equation  (3)  and  its  corresponding  boundary  equations 
(6)  and  (9)  to  find  a  discrete  solution  in  the  "inclined  advance" 
solving  region  by  transforming  the  initial  plane  no  another 
plane  (which  is  perpendicular  to  the  body  axis  of  the  rear  cone). 

Of  course,  the  initial  plane  and  the  solving  plane  must  be 
spatially  directed.  (2)  The  flow  parameters  and  shock  wave  shape 
parameters  are  converted  from  the  (z^,  r<|,  coordinate  and  the  (Z2 
’  ^2’  ^2^  coordinate  (see  Figure  1).  (3)  The  transformed  plane 

is  used  as  the  initial  plane  to  find  the  solution  at  the  object 
end  by  "parallel  advance".  The  specific  method  is  shown  in 
reference  C2].  (4)  In  the  "advance"  calculation,  we  will 

encounter  some  points  on  the  ridge  of  the  surface.  The  flow  may 
expand  or  contract  at  these  points.  Therefore,  it  is  necessary 
to  solve  these  singular  points  individually. 


Figure  1.  Definition  of  Attack  Angle  and  Sideslip  Angle 

1 .  initial  value  plane  ne 

2.  transformation  plane  n' 

3 .  shock  wave  plane 

1.  Calculation  of  Points  in  the  Flow  Field 
T  is  used  to  express  a  column  vector  whose  components  are  p, 
p,u,  V,  and  w.  For  equation  (3),  the  pre-estimation  step  and  the 
modification  step  of  the  MacCormack  scheme  can  be  written  as 


7 


•  ♦  1 


(10) 

(11) 


The  values  of  af /aC  and  af /av  ,  which  are  derivatives  of  flow 
parameters  with  respect  to  ^  and  ^  included  in  the  function  of 
af/an>  were  approached  by  the  forward  (pre-estimate)  and  backward 
(modification)  derivation  difference.  When  the  radial  mesh 
spacing  AC  and  circumferential  mesh  spacing  Atp  are  ascertained, 
the  integration  step  An  is  determined  by  the  stability  condition. 

2.  Calculation  of  Shock  Wave  Points 

The  differential  equation  (6)  was  solved  by  discretization 
using  the  MacCormack  scheme.  The  pre-estimate  and  modification 
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steps  are: 


^-411  (12) 


The  pre-estimate  and  mcdlflcatlon  steps  for  the  shock  wave  /36 
position  are: 


(i/i)::i=(/.);.  *  +  (m,);,  (14) 

(15) 

The  values  of  dm^/a«r  in  the  right  function  of  am^/dn  can  be 
approached  by  the  central  difference  of  the  values  at  various 
9,  The  values  of  af/ac  for  various  derivatives  of  flow 
parameters  with  respect  to  9  are  approached  by  the  same  method  as 
that  for  Internal  points.  The  values  of  af/a^  are  approached  by 
a  two  point  backward  difference  method  in  the  pre-estimation 
step.  In  the  modification  step,  they  are  approached  by  a  three 
point  derivation  difference  method. 

The  simple  process  of  shock  wave  point  computer  is  described 
as  follows.  The  shock  wave  slope  m^  was  determined  by  the 
difference  scheme  (12)  or  (13).  The  shock  wave  position  f^  was 
obtained  from  the  difference  scheme  (14)  or  (15).  fJj^  was 
determined  by  the  central  difference  of  the  f^  value  on  each 
meridian  plane.  Then,  the  circumferential  slope  of  the  shock 
wave  f^^  and  m2  can  be  obtained  from  the  following  equation: 


In  the  above  equation 


E,+/:,Er 


l-m,  E, 


(16) 


Finally,  flow  parameters  behind  the  wave  were  determined  by  m^ 
and  m2  using  the  Ranke-Hugonlot  relation. 

3.  Computation  of  Surface  Points 

The  same  difference  scheme  as  that  for  internal  points  is 
used  to  discretize  and  solve  the  surface  point  equation  (9).  In 
the  modification  step,  the  radial  differential  af/d?  included  in 
the  right  function  is  approached  by  a  three  point  deviation 
difference  method. 

4.  Tightening  of  Mesh  Points  Near  the  Surface 

A  boundary  surface  =  const  was  set  up  to  divide  the 
meshes  near  the  surface.  Smaller  meshes  are  used  inside  and 
larger  meshes  are  used  on  the  outside.  If  the  total  number  of 
radial  meshes  is  N  and  the  boundary  is  set  up  at  i  =  n' ,  then  the 
mesh  spring  inside  (close  to  the  surface  of  the  object)  is 

Outside  the  boundary,  it  is 


Therefore,  the  meshes  near  the  surface  are  significantly 
tightened  as  compared  to  uniform  meshes  while  the  total  number  of 
meshes  does  not  increase. 

IV.  Calculation  of  Flow  Characteristics  at  Singularity  Points 

on  the  Surface 

Because  the  flow  parameters  on  the  intersecting  line  of  the 
front  and  rear  cones  have  multiple  values,  therefore,  they  must 
be  solved  individually.  The  specific  method  is  to  establish  a 
right  angle  coordinate  ,  t2 ,  F^)  at  the  origin  A  on  the 
intersect.  F2  is  the  direction  of  the  unit  normal  vector  of  the 
front  cone  surface,  F^  is  in  the  direction  of  the  tangent  of  the 


inner  intersecting  line  and  t^  is  perpendicular  to  both  t2  and  t^. 


7  '  —  t  t  i 

* '  I  7.x  7il 

(D) 

The  supersonic  flow  expands  or  contracts  on  the  plane  /37 

so  that  the  flow  direction  is  shifted  by  a  certain  angle.  The 
flow  parameters  satisfy  the  two-dimensional  Prandtl-Meyer 
solution  or  two-dimensional  shock  wave  correlation.  Moreover, 
the  velocity  component  along  t^  does  not  change  beyond  point  A. 


V.  Numerical  Examples 

Several  numerical  examples  are  given  in  the  following. 

Figure  2  shows  the  pressure  distribution  around  a  bent  cone 
obtained  by  the  "parallel  advance"  and  "Inclined  advance"  method. 
The  "inclined  advance"  method  began  with  the  initial  value  plane 
(0,  =  90“ )  to  =  74“ ,  then  from  0^  =  74“  to  0*  =  90“ .  From 
the  figure  we  can  see  that  the  results  of  these  two  techniques 
are  identical.  The  same  accuracy  resulted. 

Figures  3  and  4  show  the  pressure  distribution  on  the 
surface  of  the  object  around  bent  cone  A.  The  geometric 
parameters  are  0^  =  10“  ,  0a  =  7“  ,  n  =8  and  0j^  =  6“ . 

Furthermore,  there  is  a  transition  link  between  the  front  and 
rear  cones.  Therefore,  there  are  two  intersecting  lines  on  the 
surface.  When  the  angle  of  attack  is  large,  the  inviscid  flow 
equation  cannot  appropriately  describe  the  true  flow  pattern. 
Furthermore,  the  attack  angle  of  the  front  cone  ®  ®N 

in  the  computation.  Therefore,  the  attack  angle  range  is  greatly 
related  to  the  degree  of  bending  of  the  bent  cone.  In  this 
example,  the  largest  angle  of  attack  was  14”  and  the  results  are 
still  good. 


Figure 


Figure 


2. 


Comparison  of  Pressure  Distributions  on  Bent  Cone 
by  "Parallel  Advance"  and  "Inclined  Advance"  Methods 

1.  parallel  advance  solution  [2] 

2.  inclined  advance  solution 


Surface  Pressure  Distribution  of  Bent  Cone  A 
1.  transition  section 
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Figure  4.  Surface  Pressure  Distribution  of  Bent  Cone  A 
1.  transition  section 

Figures  5,  6,  and  7 -are  the  flow  field  pressure 
characteristics  and  Internal  shock  wave  shape  of  bent  cone  B  (e^ 

=  e,  »  7* ,  0JJ  =  13”)  with  sideslip.  The  mesh  numbers  used  in  the 
calculation  were  M  =  20  (radial)  and  M  =  20  (circumferential). 

No  filtering  was  used  to  smooth  the  process.  The  inner  shock 
wave  position  in  the  figure  is  determined  by  the  pressure 
variation  in  the  flow  field.  From  Figure  7  we  can  see  that  the 
pressure  change  position  (i.e.,  a  strong  Interruption  position) 
is  very  obvious  when  an  inner  shock  wave  appears.  Therefore,  it 
is  basically  correct  to  determine  the  inner  shock  wave  position 
based  on  it.  In  addition,  the  result  of  total  entropy  was 
compared  with.  From  the  numerical  examples,  the  calculated 
entropy  of  flow  field  points  and  boundary  points  is  in  agreement 
with  that  of  the  incoming  flow. 


Figure  5 


Shape  of  Shock  Wave  of  Bent  Cone  B  (expansion 
wave  not  shown) 


H 

3 

i 


1 .  external  shock  wave 

2.  inner  shock  wave 


(za  13.4.> 

Figure  6.  Shape  of  Shock  Wave  of  Bent  Cone  B 

1 .  external  shock  wave 

2.  inner  shock  wave 


Figure  7,  Pressure  Characteristics  at  Z  =  14.4  for  Bent 
Cone  B 
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A  NUMERICAL  COMPUTATION  FOR  THE  INVISCID 
SUPERSONIC  FLOW  AROUND  BENT  CONE 


li  Chuqun 

(Beijing  Institute  of  Aerodynamics) 

Abstract 

A  numerical  computation  method  for  the  inviscid  supersonic  flosr 
around  a  bent  cone  is  given  in  this  paper  by  using  finite-difference  method 
and  shock-capturing  technique.  In  order  to  overcome  the  difficulty  due 
to  the  bend  of  the  bent  cone  axis  in  the  advanced  process  solving  equation, 
the  solving  plane  of  Euler  equation  is  considered  as  a  inclined  plane 
with  the  advanced  direction  using  a  transform  of  the  equation  indepent- 
dent  variables.  The  second  order,  two-step  MacCormack  scheme  is  used. 
The  equation  system  for  boundary  compution  is  derived  according  to  the 
characteristic  compatibility  relations  on  body  surface  or  shock  wave  and 
corresponding  boundary  condition.  The  flow  properties  of  the  singular 
points  on  the  intersected  line  between  two  cone  surface  are  solved  indivi¬ 
dually.  The  non-uniform  computation  meshes  along  the  direction  of 
body  radius  are  used  in  order  to  raise  the  computation  accuracy  when 
the  gradient  of  the  entropy  near  body  surface  is  very  large.  The  compu¬ 
tation  results  in  many  numerical  examples  indicate  that  good  results 
can  be  obtained  by  means  of  this  method. 
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Theory  of  Similarity  and  the  Profile  of  the  Mean  Density 
Distribution  of  Hypersonic  Boundary  Layer 
Cal  Shutang  (Department  of  Modern  Mechanics,  The 
University  of  Science  and  Technology  of  China) 

Gao  Shouen  (Department  of  Physics,  Hangzhou  Teachers  College) 

Abstract 

It  is  well  known  that  in  a  boundary  layer,  the  variation  of  [ij 
the  physical  quantity  along  the  flow  direction  (the  x  direction) 
is  much  smaller  than  that  along  the  direction  perpendicular  to 
the  flow  (the  y  direction).  At  the  same  time,  the  velocity  in 
the  "y"  direction  is  much  smaller  than  that  in  the  the  "x" 
direction.  In  this  paper,  we  derive  the  vorticity  pulsating 
equations  by  means  of  neglecting  the  physical  quantity  variation 
in  the  "x"  direction  and  the  velocity  variation  in  the  "y" 
direction.  In  a  coordinate  system  of  the  average  fluid  velocity 
motion,  we  introduce  the  similarity  hypothesis  to  determine  the 
pulsating  velocity  and  the  pulsating  density.  Based  on  the 
necessary  condition  of  the  pulsating  similarity,  we  find  the 
equation  p  =  k(y  +  yQ)”  as  the  average  density  variation  rule  in 
the  turbulent  boundary  layer  of  the  compressible  fluid.  Although 
this  method  is  very  approximate,  it  conforms  to  the  experimental 
data.  This  paper  demonstrates  that  in  a  comparatively  large 
range,  the  similarity  hypothesis  is  applicable. 

I.  Introduction 

Many  papers  have  been  published  about  the  turbulent  boundary 
layer  of  the  incompressible  fluid  on  a  flat  surface.  Although 
some  other  papers  discussed  the  turbulent  boundary  layer  of  the 
compressible  fluid  on  a  flat  surface,  yet  their  discussions  were 
incomplete,  and  their  selections  of  the  close  systems  were  quite 
arbitrary.  In  this  paper  we  adopt  the  vorticity  similarity 
conception  to  study  the  profile  of  the  mean  density  distribution 


of  the  hypersonic  boundary  layer.  We  directly  applied  Prof. 

Zhou  Peiyuan's  method  (published  in  1959)  of  the  compressible 
fluid  vorticity  equations.  Based  on  the  premise  of  the 
similarity  conditions,  we  obtain  a  few  groups  of  the  related 
physical  quantity  equations.  They  can  be  simplified  to  satisfy 
the  differential  equations  for  the  mean  density.  After 
integrating  we  obtain  an  exponent  equation  for  the  mean  density 
profile.  This  process  seems  to  be  too  approximate.  However,  the 
mathematical  results  conform  with  the  experiment  results  in  a  wide 

range. 

II.  Vorticity  Pulsating  Equations  of  Compressible  Fluid 
First,  let  us  put  down  the  Navier-Stokes  equations  for 
compressible  fluid  and  the  Reynolds  equation  after  the  Navier- 
Stokes  equations  being  averaged  by  the  Reynolds  method,  i.e. 


( 

p  (/,+ P^JU; + iron’) 


_  dp  .  a  - 

dxf  dxi  "  (2) 

where,  *  velocity,  p'  =  pressure,  p  =  density,  =  stress 
tensors  other  than  pressure.  As  to  the  relation  of  the  mean 
quantity  and  instantaneous  quantity,  we  have 


where.  A*  =  physical  quantity  and  A  =  pulsating  quantity.  When 
we  substract  (1)  from  (2),  we  obtain 

-^(p«; +£7.p'+p'".'  -p'ttJ) £7.“) +  p +p'(7'C/» 


+  P  Ui  b;  -  p  u;  U'l  ~U'P'v'i  +  0  iP’“i  ~  U  ip'u  'i 

+  P'U,'u',  —  fl'Biu'i)  “  -4^ - 


similarly,  we  have  a  continuous  equation  and  an  averaged 
continuous  equation: 


(4 

(5 

The  subtraction  between  these  two  equations  results  in  the 
pulsating  density  equation: 


In  equation  (3),  when  we  retain  the  spin  factor,  cancel  "p" 
simplify  it  with  equations  (5)  and  (6),  and  let 


fi>«* 


duj  _  duj 
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a* 


dxt  dx, 


then  we  can  obtain  the  equation  to  satisfy  the  pulsating 
vorticity 


=  d 
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in  which. 


•  T* 


III.  similarity  and  Mean  Density  of  the  Supersonic  Boundary 
Layer 

First,  let  us  consider  the  steady  state  conditions,  In  which 
the  partial  derivative  of  the  average  quantity  of  equation  (7) 
with  respect  to  time  Is  zero.  Secondly,  when  we  consider  the 


supersonic  boundary  layer,  the  "x"  derivative  term  of  an  average 


quantity  can  be  neglected  as  compared  to  Its  ”y"  derivative. 


At 


the  same  time  we  can  omit  the  velocity  terms  In  the  "y" 
direction.  Then  by  adopting  Zhou's  method,  we  Introduce  a 


coordination  system  of  motions  which  varies  with  the  point 


and  also  Introduce  a  similar  condition 

ii 


Let, 


A 


then 


(8) 


In  which  "x^",  which  can  usually  be  considered  as  the  origin.  Is 
the  coordinate  of  point  "Pq"  In  the  coordinate  system  of  motion, 
"x  "  Is  the  coordinate  of  a  certain  point  In  this  system,  and  "a" 
Is  the  characteristic  length  of  the  pulsating  quantity.  We 
further  substitute  (8)  to  the  vortlclty  pulsating  equation  (7), 


and  eliminate  the  viscosity  terms  in  (7).  In  this  paper,  we  are 
not  going  to  write  the  details  of  the  derivation  because  of  the 
complexities  of  these  equations.  Here  we  only  present  the 
coefficients  in  front  of  the  nondimensional  of  all  the 

terms  in  these  equations.  There  are  seven  different  coefficients 
as  follows: 


®  ’-iff-' 


In  order  to  obtain  the  similarity  of  the  pulsating  velocity  " 
and  the  pulsating  density  "p*”  in  the  profile  of  the  boundary 
layer  (l.e.,  in  the  profile  of  the  boundary  layer,  the  functions 
of  and  are  identical),  the  ratios  of  these  seven 

constants  have  to  be  Independent  of  the  "y"  coordinate. 
Therefore,  coefficient  (3)  is  directly  proportional  to 
coefficient  (4),  i.e.. 


=  9  dt/ 


A  dy 


ccq 


dp  dU 
dy  dy 


After  eliminating,  idU/dy)  q  we  obtain 


dy  A 


(9) 


Similarly  coefficient  (5)  is  directly  proportional  to  (6),  then 


2 

When  we  eliminate  "q  ",  we  obtain 


A  dy  *  dy* 


dp 


Acc  - 
dy* 

When  we  eliminate  "a"  in  (9)  and  (10),  we  obtain 


(10) 


in  which  "C"  is  a  constant.  Simplifying,  we  obtain 


<f*p  dfi 

AL  P  (11)’ 

dy 


After  integrating,  we  obtain 


la-^-Clnp +lnif 


(12) 


in  which  InK  is  an  integration  constant. 


After  shifting  the  terms,  we  obtain 


Simplifying , 


we  obtain 


(12)’ 


^Kdy 


Then  after  integrating,  we  obtain 

in  which  Is  an  integration  constant.  After  simplifying,  we 

obtain  ,  , 


Let 


n 


1 


we  obtain 
Further,  let 

After  dividing  with  "p^" , 


P»K'(y+y,)' 


’■-S- 


we  obtain 
-^•K-(y+y,r 


(13) 


(14) 


This  non-dimensional  equation  expresses  the  mean  density. 

IV.  Comparison:  Theory  and  Experiment 
Let  us  use  the  graphic  method  to  find  the  accuracy  of  the 
equation  by  using  the  x  abclssa  to  represent  the  data  calculated 
from  the  equation  (14)  and  the  ordinate  to  represent  the 
experimental  data.  If  the  calculated  results  agree  with  the 
experiment  results,  then  all  the  corresponding  points  should  fall 
on  a  45°  line  from  origin.  The  comparisons  are  displayed  In 
Figures  1,  2  and  3.  We  can  see  that  all  the  points  are 
practically  falling  on  the  45°  line  In  each  figure.  These 
figures  explain  that  except  on  both  ends  of  each  line  (which  show 
some  deviation)  the  theoretical  calculated  data  essentially  agree 
with  the  experimental  results. 

V.  Discussion 

It  Is  concluded  that  the  similarity  assumption  Is  quite 
approximate,  yet  it  is  practical  for  application. 


Figure  1.  Po  =  3200cmHg,  a  =  4  inch,  n  =  0.98,  no=  0.15, 
K'  =  0.606 

1.  calculated 

2.  experimental 
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The  Engineering  Numerical  Technique  for  the  Determination 
of  the  Inviscid  Flow  Field  Heating  Rate  on  Ballistic 

Re-entry  Vehicles 
Yang  Maozhao,  He  Fangshang 
(China  Aerodynamic  Research  and  Development  Center) 

Abstract 

RVSPHR  (The  Ballistic  Re-entry  Vehicle's  Shock  Pressure  and 
Heating  Rate  Computer  Code)  is  an  approximate  numerical  technique 
for  the  determination  of  the  inviscid  aerodynamics  and  heat 
transfer  of  re-entry  vehicles.  In  the  transonic  flow  region 
surrounding  the  nose  stagnation  point,  the  distribution  of  the 
surface  pressure  and  the  shock  angles  is  assigned  from  the 
correlation.  With  the  assumptions  for  the  profiles  of  pressure 
and  normal  velocity  along  body-normals  in  Von  Mises  coordinates, 
local  iteration  is  carried  out  at  each  normal  to  determine  the 
shape  of  the  shock  wave  based  on  mass  continuity  satisfied  along 
that  normal.  In  the  downstream  supersonic  regions,  the  exact 
Euler  equations  of  motion  are  integrated  using  the  finite- 
difference  method.  The  surface  heating  rate  is  calculated  using 
a  non-inviscid  flow.  All  calculated  results  are  compared  with 
more  rigorous  solutions  as  well  as  experimental  results  and  good 
agreements  have  been  found. 

I.  Introduction 

The  determination  of  the  aerodynamic  characteristics  using  the 
numerical  technique  of  inviscid  flow  and  the  analysis  of  the 
aerodynamic  heating  process  using  numerical  calculation  of  the 
boundary  layer  have  been  successfully  carried  out.  However, 
since  the  rigorous  numerical  calculation  is  time  consuming  and 
costly,  it  has  not  been  widely  employed  for  carrying  out  the 
initial  engineering  estimations  in  certain  applications.  For 
engineering  application,  we  have  used  engineering  approximate 
solutions  for  treating  the  transonic  region  in  rigorous  inviscid 
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flow  calculation.  The  iterative  numerical  procedure  in  solving 
the  elliptic  equations  for  the  subsonic  region  surrounding  the 
nose  stagnant  point  is  avoided.  This  simplifies  greatly  the 
problem  of  the  subsonic  region  and  accelerates  the  calculation 
without  sacrificing  much  of  the  accuracy  or  affecting  the 
solution  of  the  downstream  supersonic  region.  With  this 
approach,  we  have  determined  the  correlation  equation  of  the 
surface  pressure  distribution  and  the  shock  angles  in  the 
transonic  region  around  the  nose  as  well  as  the  profiles  of 
pressure  and  normal  velocity  along  body-normals  in  the  shock  wave 
layer.  The  position  of  the  shock  wave  corresponding  to  each 
normal  on  the  surface  is  determined  by  local  iteration  or  overall 
iteration  based  on  mass  conservation.  The  solution  of  the 
transonic  region  provides  the  preliminary  data  of  the  profiles 
for  calculating  the  downstream  supersonic  region.  The  exact 
Euler  equations  of  motion  for  the  supersonic  region  are 
integrated  using  a  numerical  technique  with  the  finite-difference 
method. 

Although  this  technique  for  determining  the  inviscid  flow 
field  is  an  approximate  method,  it  is  a  technique  for  calculating 
the  overall  flow  field  in  regarding  the  estimation  of  pressure 
distribution  and  the  shape  of  the  shock  on  the  aircraft.  It 
provides  the  necessary  data  for  boundary  layer  calculations  with 
consideration  of  the  effects  of  entropy  swallowing.  Based  on 
this  technique,  we  have  further  derived  a  method  of  integration 
with  the  momentum  in  the  boundary  layer  and  a  method  of  entropy 
swallowing  calculation  relating  to  mass  balance  of  the  flow.  The 
new  technique  does  not  need  to  carry  out  iterations  on  the 
parameters  of  the  outside  edge  of  the  boundary  layer  and  much 
less  computer  time  is  required  in  calculating  the  boundary 
parameters  and  the  surface  heating  rate. 

We  have  carried  out  calculations  on  many  aircrafts  and  wind 
tunnel  test  models.  The  results  are  in  good  agreement  with 
these  obtained  from  more  rigorous  solutions  and  experimental 
data.  There  are  also  significant  savings  in  computer  time. 
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II.  Summary  of  the  Technique 
(1)  Technique  for  the  Determination  of  the  Inviscid  Flow 

The  purpose  of  this  technique  is  to  avoid  the  iterative 
calculation  in  solving  the  elliptic  equation  in  the  subsonic  flow 
region  surrounding  the  front  edge  of  the  passive  nose.  The 
following  approximations  and  assumptions  are  made  in  the 
calculation  of  the  transonic  flow  region  around  the  nose: 

1.  The  pressure  distribution  in  the  transonic  region  (from 


stagnation  point  to  the  sonic  point)  is  determined  using 

[1] 


engineering  equations.  Love's  modified  equation  is  used  in 
the  calculation  of  the  upstream  region  of  the  sonic  point. 

Matched  Newtonian  equations  with  Prandtl-Meyer  expansion  are  used 
for  the  region  after  the  sonic  point.  The  equations  by  curve 
fitting  for  the  latter  are  shown  in  the  following: 

For  M  <1.2, 


/»-/),{l  +  1.08sin’C(jr/2)-fl)3[exp(-4>/  In  s*/s  ) -l]-0.468  exp(-4v'ln  i*/s)} 

2.1-1 

For  M  >1.2, 

2.1-2 

Where  0,  s,  0  ,  s  are  the  dip  angle  and  arc  length  of  the 
surface  at  local  point  and  the  sonic  point  respectively.  Pq  is 
the  pressure  at  the  stagnation  point,  which  can  be  expressed  by 
the  following  equation  with  ideal  gas: 


{y±L  Ml  y 

'2  /  \2vMi-y  +  l/ 


2.1-3 


For  equilibrated  air,  it  should  be  calculated  by  iteration 
from  the  relations  of  the  thermodynamic  characteristics  using  the 
entropy  behind  the  normal  shock  wave  and  the  accurate  enthalpy,  h 
,  at  the  stagnation  point.  All  quantities  used  in  this  paper  are 
dimensional  and  in  metric  system. 


2.  In  the  shock  wave  layer  with  relatively  small  spherical 
angle,  the  profile  of  the  normal  velocity  on  the  surface  other 
than  the  flow  at  the  stagnation  point  can  be  correlated  with  the 
profile  at  the  stagnation  point  using  the  following  equation: 


v,U) 

Vc(y.) 


2.1-4 


where  v, (5)  is  the  normal  velocity  distribution  of  the  flow  at 

the  stagnation  point,  v,  (y  ),  v((j,  )  are  the  normal  velocities 

s  s 

behind  the  shock  wave  at  and  away  from  the  stagnation  point 
respectively.  The  values  are  fixed  with  given  shape  of  the  shock 
wave.  Vj (5)  can  be  determined  using  the  pressure  distribution  of 
the  flow  at  the  stagnation  point.  According  to  the  numerical 
technique  for  the  determination  of  the  inviscid  flow  ,  Po (c) 
can  be  fitted  in  the  following  equation: 


2.1-5 


where  is  the  pressure  behind  the  normal  shock  wave.  From 

Equation  (2.1-5),  the  enthalpy  distribution,  hjC^),  can  be 
calculated  and  the  following  equation  is  obtained: 


gJZh.-h,U)^ 


2.1-6 


In  Von  Mises  coordinates,  the  relationship  of  the  flow 
function  and  the  normal  coordinate  y  is 


urdy 


which  gives 


2.1-7 


2.1-8 


In  the  region  near  the  stagnation  point,  using 

f  ^yJy,.i  =  yly,  2  19 

the  5  value  corresponding  to  each  point  y  in  the  flow  at  the 


stagnation  point  can  be  easily  calculated. 

3.  The  following  profile  of  pressure  along  body-normals  in 
the  shock  wave  layer  is  assumed: 

=  p»+ 2  o  1-10 

where  p^  is  the  estimated  surface  pressure  from  previous 
engineering  equations,  p^j^  is  the  value  of  the  pressure  behind 
the  shock  wave  in  the  profile. 

4.  The  initial  shape  of  the  shock  wave  is 

!/«»  .=yi»  ,.,dx  2.1-11 


2.1-12 


=  0.5326 -i- 0.3333  +  0.2122  fli.i 


2.1-13 


The  above  calculations  can  be  carried  out  if  the  distance,  ygjj»o 
from  the  body  to  the  shock  wave  on  the  flow  line  at  the 
stagnation  point  is  available.  The  value  of  ygj^»Q  can  be 
calculated  from: 


2.1-14 


For  an  ideal  gas  the  following  equation  may  be  used 


2.1-15 


After  obtaining  the  values  of  the  flow  field  parameters,  p  and  u 
of  the  calculating  cross-section,  the  new  position  of  the  shock 
wave  y^  can  be  determined  by  the  requirement  of  mass  continuity 
which  gives 


2.1-16 


Since  the  calculated  position  of  the  shock  wave  may  not  be 
consistent  with  the  original  positicm  assigned ,  additional 
iterations  are  necessary  until  the  fl|lults  of  consecutive 
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iterations  show  the  desired  accuracy.  This  is  the  local 
iteration  process  for  obtaining  the  shape  of  the  shock  wave. 

With  this  calculation  technique,  we  have  also  developed  an 
overall  iterative  procedure.  After  obtaining  the  positions  of 
the  shock  wave  on  the  cross-sections  along  all  body-normals  by 
local  iterative  procedure,  the  new  distribution  of  the  shock 
angles,  can  be  obtained  by  parabolic  fit  or  spline  fit  of 

the  shape  of  the  shock  wave.  Equation: 

is  used  to  replace  the  old  equation  of  the  shock  angles  (2.1-13) 
and  a  new  local  Iteration  for  the  profile  along  each  body-normal 
is  carried  out  for  further  improving  the  estimation  of  the 
position  of  the  shock  wave.  Better  convergency  of  the 
calculation  can  be  obtained  with  62  =  0.2~0.6. 

The  above  calculation  is  carried  out  from  one  cross-section 
to  next  cross-section  starting  from  the  point  near  the  stagnation 
point  until  the  velocity  u  reaches  the  sonic  velocity  at  one 
particular  cross  section.  After  that  the  calculation  should  be 
carried  out  following  the  procedure  for  supersonic  region. 

Since  the  calculation  in  the  supersonic  region  is  carried 
out  in  boundary  layer  coordinate  system  (s,y),  which  is  different 
from  the  Von  Mlses  coordinate  system  (  s^,4>)  and  the  preliminary 
cross-section  data  for  the  calculation  should  be  obtained  by 
interpolation,  the  simultaneous  equations,  as  shown  in  the 
following  vectors,  for  calculation  of  the  supersonic  region  can 
be  derived  from  Euler  equations  of  motion  for  the  stable, 
inviscid  and  compressible  axisymmetric  body. 

^  +  2.1-18  /59 


where 


A—  1  /Ao  +  Ca  H  +  ari  dy,  \  A  +  C 
pu^  \  £“  y,  dx  )  E 


(sin  6  +  a  cos  0) 


£^\sintf^acos5)-i±^ 


£ 

if-H-Q///?) 


a^rj  dy, 
u^y,  dx 


>»  ^vly.ix) 


E  =  l-(a‘lu‘) 


o  =v/u 


r.r.+ycose 


At  point  C  =  (n  +0)a?,  n  =  Cm  +  (1/2)]An,  e  =  0.625  the  above 
equations  can  be  rearranged  in  implicit  forms  to  give: 


where 


r Cm*l  +  A  £-=77 

***7 


2.1-19 


/7...=[i-2<i-  0  J-2(i-  e  >  2ak:| 


and  I  is  the  unit  matrix  .  iVro  of  the  elements  can  be  expressed  as 


1 


2yJ 


2*;;]=  $  u::i +«:**) +a-  «)<«:♦,+*:) 


The  solution  of  the  above  finite-difference  equations  can  be 
obtained  by  the  catching  up  technique.  The  recurrence  relations 
of  the  coefficients  are  as  follows: 


2.1-20 


2.1-21 


where  is  the  inverse  of  ,  I  I- -  I  I  represents  the 

determinant. 

In  addition 


With  impermeable  body  surface,  o=  0,  then 


2.1-22 


where 


2.1-23 


(7,»o 


Therefore,  the  values  of  ,  Go  within  the  shock  wave  can  be 
determined  by  the  recurrence  relations  starting  from  the  surface 
of  the  body.  Combining  with  the  shock  wave  condition  of  Rankine' 
Hugoniot,  is  readily  obtained  from  the  relation  of  the  shock 


Then,  the  values  of 

Cif-i,  Car. I,  . ,  if 

can  be  obtained  from  the  finite-difference  equation  (2.1-19)  by 
calculating  back  to  the  surface. 


II.  Improved  Entropy  Swallowing  Technique 
It  is  well  known  that  the  effects  of  entropy  swallowing  of 
the  boundary  layer  are  usually  determined  by  mass  balance.  The 
mass  flow  swallowed  by  the  boundary  layer  should  be  equal  to  the 
mass  flow  which  passes  through  the  bow  shape  shock  wave  along  the 
inviscid  flow.  The  equation  of  mass  conservation  of  the  flow 
path  can  be  written  as 

M-  r‘,  =p,  II,  r 

^  2.2-1 

where  the  momentum  thickness  of  the  boundary  layer,  e., ,  can  be 
determined  from  the  following  axisymmetric  equation  of  momentum: 


-JT 


'2+/f,  du,  ^  1  dp,  ^  1  dr  \ 
«.  ds  p,  ds  r  ~Ss ) 


2.2-2 


Iterations  are  needed  to  converge  the  values  on  both  sides  of 
equation  (2.2-1).  The  estimated  shock  angle  (or  shock  height  r_ ) 
as  well  as  the  values  of  the  parameters  p- ,  u.  and  M.  should  be 
modified  continuously.  The  convergency  can  be  tested  by 
substituting  the  value  of  e  determined  from  equation  (2.2-2)  in 
equation  (2.2-1).  This  is  the  iterative  process  used  previously 
for  entropy  change.  In  our  analysis,  the  axisymmetric 
integration  equation  of  momentum  is  written  in  terms  of  the  total 
differential  form  of  the  mass  flow.  Equation  (2.2-2)  is 
rewritten  as 

^  (p.  II.  r  fl*)  --L  p.  II,  rc,-  (l+H.)  (p.  «.  2.2-3 


where  is  the  mode  parameter  depending  on  the  flow.  It  can  be  /61 
written  in  the  following  two  expressions  for  laminar  and 
turbulent  flows  respectively. 

For  laminar  flow: 


For  turbulent  flow: 


59156  hJh.+h.lh.-\ 
/f,«0.14286l,(l-7l,) 


I.- 

I.- 


!; 

i: 


_ _ 

A»  +  (A,  — — 

_ _ 

h.+  {h.-h,)g+ 


dz 


dz 


2.2-4 


2.2-5 
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where 

9P.u.h.  (2.2-6) 


c  can  be  related  to  the  heat  flux,  q  ,  by  Reynold’s  analogy.  It 
f  w 

should  be  noted  that  equation  (2.2-3)  is  different  from  equation 

(2.2-2).  The  mass  flow  PcVic^^M  is  considered  to  be  a  single 

parameter  and  it  can  be  calculated  by  Euler's  method  or  Runge- 

kutta  method  of  integration  utilizing  the  solution  of  the 

inviscid  flow,  tg  is  then  obtained  from  equation  (2.2-1)  and  the 

entropy  behind  the  shock  wave  at  position  tg  can  be  calculated  by 

interpolation  from  the  shape  of  the  inviscid  shock  wave.  The 

peripheric  parameters  at  the  outer  edge  of  the  boundary  layer,  p^^ 

,  u^ ,  h^ ,  can  be  calculated  along  the  flow  using  the  pressure 

distribution  of  the  inviscid  flow.  Iterative  modifications  of 

the  value  of  du^./ds  in  equation  (2.2-3)  are  carried  out.  The 

above  calculation  is  made  repeatedly  with  newly  calculated  Uj. 

value  until  the  convergency  is  reached.  Practically,  the 

convergency  can  be  reached  shortly.  It  takes  much  less  computer 

time  than  the  technique  of  variable  entropy  with  iterative 

calculation  on  all  the  peripheric  parameters  such  as  p^ ,  u^  and 

h^.Thls  is  an  improved  numerical  technique  of  variable  entropy. 

III.  Results  of  Calculation  and  Comparison 
All  calculated  results  of  the  inviscid  flow.  Including  shock 
shape,  surface  pressure  distribution  and  the  flow  field  profiles 
through  the  shock  la^er,  are  compared  with  more  rigorous 
numerical  solutions^  as  shown  in  Figures  1,  2  and  3. 

Calculations  are  carried  out  on  a  semiconic  body  having  a  cone 
angle  of  10“  at  various  Mach  numbers  and  elevations  with 
equilibrated  air  and  ideal  gas.  The  results  Indicate  that  the 
shock  shape  and  surface  pressure  distribution  are  quite 
consistent  with  vigorous  numerical  solutions.  The  position  of 
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shock  wave  determined  by  local  iteration  and  the  off-body 
distance  of  the  shock  wave  at  the  stagnation  point  obtained  by 
parabolic  interpolation  are  compared  with  rigorous  numerical 
solutions  and  show  relative  deviations  within  1.1%  to  1.8%  as 
listed  in  Table  1.  The  density  and  velocity  profiles  of  the  flow 
in  the  shock  layer  at  =  2.0  are  in  good  agreement  with 
reference  [2].  Some  deviation  of  the  pressure  profile  is 
observed  near  the  body  surface,  while  good  agreement  is 
approached  near  the  shock  wave. 


1.  our  result 

2.  numerical  solutions  of  [2] 


Figure  2(a).  Surface  Pressure  Distribution  on  the  Nose 

1 .  ideal  gas 

2.  our  result 

3.  numerical  solutions  of  [2] 
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Figure  2(b).  Surface  Pressure  Distribution  on  the  Body 

1 .  ideal  gas 

2.  our  result 

3.  numerical  solutions  of  [2] 


1.  equilibrated  air 

2.  our  result 

3.  numerical  solutions  of  [2] 


Obviously,  the  approximation  technique  for  the  transonic  flow 
region  does  not  cause  significant  error  in  overall  calculation  of 
the  inviscid  flow  field  and  much  less  computer  time  is  required 
than  more  rigorous  solutions.  This  technique  is  very  effective 
and  feasible  in  engineering  applications. 


Table  1 


_i 

T-r 

1  X«[2] 

•.05$ 

1  0.054 

0.1315 

1  0.1300 

Key  on  following  page. 
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1.  equilibrated  air 

2.  ideal  gas 

3.  our  result 

4.  reference  [2] 

We  have  also  made  a  calculation  based  on  the  experimental 
model  of  Widhopf^^^.  The  pressure  distribution  is  shown  in 
Figure  4.  The  calculated  results  are  in  good  agreement  with  the 
experimental  results.  The  calculated  heat  flux  at  various  flow 
conditions  is  shown  in  Figures  5,  6  and  7.  The  experimental 
results  of  reference  [3]  and  the  boundary  layer  numerical 
solutions  of  reference  [4]  are  also  shown  in  the  figures.  Our 
calculated  results  are  quite  satisfactory  in  comparison. 


Figure  4.  Comparison  of  the  Calculated  Pressure  Distribution 
With  Experimental  Data 

1,  our  results 

2.  experimental  data  of  C3] 
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Figure  5.  Distribution  of  Heat  Flux 

1.  our  results 

2.  numerical  solutions  of  [4] 

3.  experimental  data  of  [3] 
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Figure  6.  Distribution  of  Heat  Flux 

1.  our  results 

2.  numerical  solutions  of  [4] 

3.  experimental  data  of  [3] 
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Figure  7.  Distribution  of  Heat  Flux 

1.  our  results 

2.  numerical  solutions  of  [4] 

3.  experimental  data  of  [3] 

Our  numerical  technique  can  be  extended  to  three  dimensional 
calculation  with  angle  of  attack  using  suitable  coordinate 
transformation  and  analogical  approximation  of  axisymmetry. 
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Visualization  of  Three  Dimensional  Separated  Flows  Using 

the  Smoke  Wire  Technique 
Xia  Xuejian 

(Beijing  Institute  of  Aeronautics  and  Astronautics) 

Abstract 

The  smoke  wire  technique  is  used  for  visualizing  the  vortex 
flow  at  the  wake  of  an  axisymraetric  body  and  the  asymmetric 
vortices  which  emanated  from  the  nose  of  the  aircraft  at  high 
angle  of  attack.  The  smoke  wire  technique  can  be  employed  for 
automatic  photograpiis,  by  the  sequential  circuit,  conditional 
photographs  or  manual  operation. 

I.  Introduction 

Various  techniques  for  visualization  of  the  flows  have  been 
applied  in  aerodyamic  research,  especially  in  some  complicated 
flows  such  as  the  separated  flow,  jet  flow  and  wake  flow  emanating 
from  wings  and  the  fuselage  at  high  angles  of  attack.  They  are 
also  effective  tools  for  studying  the  separated  flow  and  the 
structure  of  the  turbulence  on  nonflying  passive  bodies. 

Recently,  various  new  techniques  for  visualizing  the  flows  have 
been  explored. 

The  smoke  wire  technique  has  been  applied  since  the  fifties. 
Recently,  this  technique  has  been  gradually  used  in  the 
visualization  of  complicated  flows.  Reference  [1]  employed  this 
technique  to  visualize  the  bubble  separation  of  the  laminar  flow 
and  the  transition  of  the  boundary  layer  of  the  flow  about  the 
wings.  It  indicates  that  the  technique  can  be  used  in  flows 
which  cannot  be  measured  by  hot  wire  probe  or  pressure 
transducer.  We  have  used  the  smoke  wire  technique  in  studying 
the  three  dimensional  separated  flows  such  as  the  waxe  flow  on  a 
three  dimensional  passive  body.  The  flow  pattern  exhibited  was 
consistent  with  the  measurements  obtained  with  hot  wires. 

Various  periodic  characteristics  of  the  flow  were  observed 


visually.  This  technique  has  also  successfully  detected  the 
presence  of  asymmetric  vortices  emanating  from  the  nose  of  the 
body  at  high  angles  of  attack  without  lateral  gliding  in  wind 
tunnel  experiments.  Since  the  smoke  lines  are  very  thin,  detail 
flow  pattern  can  be  observed.  This  technique  can  be  widely  used 
in  research  for  its  convenience  and  less  environmental  pollution. 

The  limitation  of  this  technique  is  on  the  wind  speed  of  the 
test.  It  cannot  be  applied  at  high  Reynolds  Number.  Reference 
[2]  has  compared  in  detail  the  characteristics  of  the  separated 
flows  on  wings,  fuselage  and  whole  body  at  high  angles  of  attack 
tested  in  water  tunnel,  wind  tunnel  and  test  flight  at 
corresponding  limits  of  the  Reynolds  Number.  It  also  analyzed 
the  effect  of  Reynolds  Number  on  the  separated  flows.  Since  the 
separated  flows  around  sharp  edges  or  fixed  edges  are  less 
sensitive  to  Reynolds  Number,  the  smoke  wire  technique  can 
provide  more  quantitative  measurement  on  the  flow  pattern. 
Additional  analyses  are  necessary  for  more  complicated  flows 
which  are  affected  by  the  Reynolds  Number.  The  flow  spectrum  one 
fuselage  at  high  angles  of  attack  using  smoke  wire  technique  can 
provide  qualitative  information. 

II.  The  Equipment  for  Smoke  Wire  Display  and  the  Technique  of 

Operation 

Figure  1  shows  the  equipment  of  smoke  wire  display  for 
studying  the  wake  flow.  The  system  consists  of  three  positions: 


the  first  portion  is  the  smoke  generator,  which  is  located  in  the 
wind  tunnel.  The  second  portion  is  the  process  electrical 
system,  which  includes  power  source  for  heating  the  smoke  wire 
and  the  synchro  control  box.  Thc'.  third  portion  is  the  light 
source,  camera  and  the  activation  system  (the  portion  in  the  dark 
shown  in  the  picture). 
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Figure  1 .  Equipment  for  Smoke  Wire  Display 

1 .  The  Smoke  Wire  Generator  and  the  Operation 

Figure  2  shows  the  smoke  wire  generator  which  is  composed  of 
resistance  wire,  supporting  frame,  the  feedback  passage  of  the 
resistance  wire.  The  resistance  wire  is  for  generating  the  smoke 
by  heating.  It  can  be  a  tungsten  wire,  nichrome  wire  or 
stainless  steel  wire.  Stainless  steel  wire  is  more  favorable 
since  it  has  high  resistivity,  fixed  strength  and  less  tendency 
to  become  brittle.  The  diameter  of  the  resistance  wire  depends 
on  the  experiment.  A  size  of  0.15mm  is  suitable  in  general.  For 
certain  flow  conditions  which  may  be  distributed  by  the  wire,  the 
selection  of  the  diameter  should  be  more  careful.  Reference  Cl] 
provides  some  information  on  the  selection  of  the  wire.  It 
should  be  noticed  that  the  diameter  of  the  resistance  wire 
affects  the  smoking  time.  Since  the  diameters  of  the  oil  drops 
attached  on  a  thin  wire  are  smaller,  the  duration  of  smoking  will 
also  be  shorter.  The  tension  on  the  resistance  wire  will  also 
affect  the  mass  flow  of  the  smoke.  The  equipment  employs  a 
tension  bolt  to  adjust  the  length  of  the  resistance  wire. 


o 

Figure  2.  The  Structure  of  the  Smoke  Generator 

1.  resistance  wire 

2.  top  bracket 

3.  bottom  bracket 

4.  top  support 

5.  bottom  support 

6.  supporting  pole 

7.  tension  bolt  for  the  steel  wire 

8.  tension  bolt  for  the  steel  wire 

9.  receptacle  for  conducting  ’*ire 

Specially  selected  oils  should  be  used  for  smoke  generation. 
Smoke  oil  for  model  trains,  paraffin  oils  as  well  as  other  oils 
have  been  used  in  experiments.  Experimental  results  indicate 
that  glycerine  can  provide  longer  smoking  duration  and  more  clear 
flow  pattern.  The  smoking  duration  can  be  a  few  seconds  with  a 
wind  speed  of  about  2  M/sec.  The  duration  is  long  enough  for 
taking  photographs. 

Various  tests  have  been  carried  out  for  selecting  a  better 
method  to  coat  the  oil.  The  oil  can  be  extruded  from  the  top 
bracket  by  compressed  air  or  dripped  down  from  the  top  end  of  the 
resistance  wire.  It  is  found  by  experience  that  uniform 
distribution  of  oil,  which  gives  uniform  flow  spectrum,  can  be 
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obtained  by  coating  manually  on  the  resistance  wire  with  a  brush 


2.  Electrical  System  of  the  Smoke  Wire  Process 


1.  flash  light 

2.  flash  light 

3.  conditional  output 

4.  auxiliary  output 

5.  smoke  wire  input 

6.  smoke  wire  output 

7.  auxiliary  output 

8.  camera 

9.  camera 

10.  relay 

11.  relay 

12.  gate  circuit 

13.  time  delay  unit 

14.  time  delay  unit 

15.  conditional  input 


16.  initiate  switch 

17.  potentiometer  -1 

18.  potentiometer  -2 

19.  rectifying  bridge 

20.  voltage  stabilizer 

21.  bypass  switch 

The  sequence  of  activating  the  smoke  wire  and  the  camera  is 
controlled  by  two  time  delay  circuits.  The  first  time  delay  unit 
controls  the  duration  of  the  electric  current  which  passes 
through  the  resistance  wire  for  heating.  The  second  time  delay 
unit  controls  the  delay  time  before  activating  the  camera.  The 
time  duration  of  both  units  is  adjusted  by  the  potentiometer. 

The  arrangement  can  be  changed  when  using  conditional  input 
(the  flow  spectrum  is  obtained  by  activating  the  flash  light  at 
the  occurrence  of  a  special  event).  Figures  3  and  4  show  the 
circuit  diagram  of  the  synchro  box  and  the  distribution  of  the 
process  durations. 


Figure  4.  Diagram  of  the  Process  Durations 

1.  automatic  signal 

2.  delay  -1 


3.  delay  -2 

4.  conditional  Input 

5.  output  of  the  flash  light  and  the  camera 

6.  auxiliary  output 

7.  position  1 

8.  position  2 

9.  delayed  duration 

10.  heating  the  smoke  wire  by  electricity 

11.  "smoke  visualized" 

12.  input-output 

13.  switch  off 

14.  switch  off 

15.  logic  1 

16.  switch  on 

17.  switch  off 

18.  Input^output 

19.  switch  off 

20.  input=output 

21.  switch  off 

There  are  three  types  of  "input"  and  "output"  for  the 
circuitry: 

1)  Power  Input  and  Output  for  the  Smoke  Wire 

The  timing  circuit  of  the  smoke  wire  is  on  when  the 
direct  current  with  variable  voltage  is  connected  to  the  smoke 
wire  input  of  the  synchro  box  and  the  output  is  connected  to  the 
smoke  generator.  The  "input"  and  the  "output"  is  connected  /74 

Internally  when  the  "bypass"  switch  is  on,  which  bypasses  the 
timing  circuit. 

2)  Conditional  Input  and  Output 

This  method  can  be  used  to  photograph  the  flow  pattern 
of  a  special  event.  The  input  signal  in  this  mode  should  be 
compatible  with  the  TTL  circuit  (integrated  circuit  of 
transistor-transistor  logic).  During  working  condition,  the 
input  voltage  should  be  low  (Figure  4)  at  the  beginning  and  jumps 
to  a  high  value  when  the  event  occurs.  This  causes  the  second 
timing  circuit  in  high  voltage  at  the  same  moment  and  provides 
output  to  activate  the  flash  light  and  the  camera.  A  low  current 
(few  mAs)  triggering  circuit  can  be  used  for  activating  the  flash 
light.  A  momentary  high  current  (up  to  4  amps,  60  volts) 
triggering  circuit  should  be  employed  to  trigger  the  camera. 
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It  should  be  noticed  that  the  conditional  photograph  Is 
usually  carried  out  In  the  dark.  The  shutter  of  the  camera  Is 
open  when  the  "Initiate  switch"  Is  on.  The  flash  light  Is 
triggered  and  the  film  In  the  camera  Is  exposed  for  recording 
when  the  conditional  event  occurs  and  the  "smoke"  Is  visible. 

The  shutter  of  the  camera  Is  closed  when  the  time  Is  out. 

Usually  the  signals  of  the  hot  wire  can  be  recorded 
simultaneously  when  the  flash  light  is  triggered.  Therefore  the 
flow  pattern  observed  can  be  compared  with  the  measurements  of 
the  hot  wire  for  studying  the  flow  structure. 

3)  Auxiliary  Input  and  Output 

This  Is  a  universal  Input  which  makes  the  system  more 
flexible  In  application.  It  can  provide  an  adjustable  auxiliary 
output  with  time  delay  for  any  signal  Input.  Usually  the  best 
flow  spectrum  can  be  observed  a  few  moments  after  the  resistance 
wire  is  heated.  The  auxiliary  Input  and  output  can  be  used  In 
such  cases  for  photographing. 

One  of  the  above  three  types  of  Input  and  output  can  be 
selected  as  needed  for  processing. 

The  principle  and  the  working  procedure  of  the  synchro 
box  of  the  process  circuit  have  been  discussed  above.  In 
addition  to  the  synchro  box,  It  needs  a  power  supply  with 
variable  voltage  for  heating  the  resistance  wire.  The  current 
should  be  about  2-3  amps  and  the  variable  voltage  should  be  In 
the  range  of  0-100  volts.  Figure  5  shows  the  circuit  diagram  of 
the  power  source  for  heating. 


Figure  5.  Circuit  Diagram  of  the  Power  Source  for  Heating 
the  Smoke  Wire 

1.  momentary  switch 

2.  switch 


3.  Light  Source  and  the  Photographic  System 

Strobe  light  was  used  as  the  light  source  on  top  of  the  wind 
tunnel  for  studying  the  wake  flow.  The  frequency  of  the  flash 
was  adjustable  from  1-13/sec.  A  plane  light  source  was  employed 
for  obtaining  the  flow  pattern  of  a  cross-section  of  the  three 
dimensional  flow.  The  camera  was  a  regular  one  of  type  135.  The 
shutter  was  activated  by  a  compressed  air  switch  controlled  by 
the  triggering  circuit  with  the  application  of  the  auxiliary 
input/output . 

Described  above  is  the  smoke  wire  technique  which  we  have 
used  to  Investigate  the  wake  flow  of  a  three  dimensional  passive 
body.  Recently,  we  visualized  the  separated  flow  with  positive 
results  at  high  angle  of  attack  using  a  simple  manual  operation. 

The  equipment  included  a  smoke  generator,  a  coupling  transformer 
and  a  resistance  wire  heated  by  AC.  Pictures  were  taken  at  the 
best  conditions  depended  on  visual  observation  with  naked  eye. 

The  smoke  wire  technique  with  manual  operation  is  simple  and  easy 
when  the  conditional  photograph  is  not  necessary  and  manpower  is 
available . 

III.  Examples  of  Visualization  of  the  Three  Dimensional  /75 

Separated  Flows  Using  Smoke  Wire  Technique 

Hot  wire  measurement  of  the  wake  flow  of  an  axisymmetrlc 
body  with  declined  base  surface  indicates  that  a  strong  periodic 
flow  appears  near  the  wake  at  certain  angles  of  decline  (such  as 
70",  60").  The  periodic  phenomenon  becomes  less  regular  with  a 
declining  angle  of  90"  (vertical  tail  cut).  The  result  has  been 
demonstrated  by  visualizing  with  smoke  wire.  Pictures  3  and  2  of 
the  flow  patterns  obtained  by  the  smoke  technique  (with  base  angle 
of  70®  and  90")  are  shown  in  the  inside  back  cover.  The  front 
end  of  the  model  was  a  semi-ellipsoid  and  the  crabbing  angle  and 
angle  of  attack  were  zero. 

Picture  1  indicates  clearly  the  regular  separation  of  the 
vortex  which  is  more  serious  at  the  top  than  at  the  bottom.  The 
same  characteristic  was  also  detected  by  hot  wire  measurement. 
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The  separation  of  the  vortex  with  a  base  angle  of  90°  as  shown  in 
Picture  2  was  less  regular  than  that  at  an  angle  of  70°,  while 
the  Intensity  was  about  the  same  from  top  to  bottom.  Intensive 
vortices  occurred  only  occasionally.  This  characteristic  is 
consistent  with  that  obtained  by  hot  wire  measurement.  The 
detailed  results  were  reported  in  reference  [3]. 

Picture  3  (Inside  back  cover)  is  the  photograph  of  the 
vortex  flow  emanating  from  the  nose  of  the  body  at  high  angle  of 
attack  obtained  by  smoke  wire  technique.  The  flow  pattern  is 
quite  similar  to  that  obtained  in  water  tunnel  testsas  reported 
in  reference  [2]  (Figure  31  in  reference  [2]). 

Figure  6  is  a  sketch  of  Picture  3.  The  model  is  at  zero 
crabbing  angle  and  with  an  angle  of  attack  of  46°.  The  picture 
shows  clearly  two  asymmetric  vortices  emanat  ing  from  the  front 
edge  of  the  nose.  The  center  portions  of  the  vortices  are  shown 
in  the  picture.  The  vortex  at  the  left  side  is  at  higher 
position  than  the  one  at  the  right  side.  The  left  vortex 
proceeded  a  longer  distance  before  breaking  apart  while  the  one 
at  the  right  side  broke  into  spiral  eddies  moving  in  reversed 
direction.  According  to  the  analyses  for  slim  fuselage  in 
reference  C2],  the  angle  to  initiate  the  asymmetric  vortices  is 
slightly  retarded  and  the  lateral  position  of  the  vortices  is 
moved  outward  somewhat  at  low  Reynolds  Number.  However,  it  has 
qualitatively  recognized  the  asymmetry  and  the  breakdown  of  the 
vortices . 


46* 

Figure  6.  Sketch  of  the  Breakdown  of  the  Asymmetric  Vortices 

at  the  Ncse  of  the  Fuselage 

1.  left  vortex 

2.  right  vortex 

3.  asymmetry  of  the  vortex 

4.  breakdown  of  the  vortex 

Picture  4  (inside  back  cover)  (from  reference  [4])  shows  the 
back  vortex  on  the  right  side  of  the  body  with  an  angle  of  attack 
of  31®.  The  swirl  and  the  position  of  the  center  of  the  vortex 
are  clearly  observed. 

IV.  Conclusion 

The  smoke  wire  technique  can  be  used  to  visualize  the  two 
and  three  dimensional  stable  flows  as  well  as  three  dimensional 
separated  flows.  It  can  provide  quantitative  as  well  as 
qualitative  information  on  the  characteristics  of  the  separated 
flows.  It  can  be  employed  to  take  conditional  photographs  of  the 
flow  pattern  at  certain  special  occasions.  The  method  is  simple 
and  flexible.  It  can  provide  the  detailed  structures  of  the 
flows.  The  smoke  wire  technique  is  an  effective  tool  for 
studying  complicated  flows. 

Picture  3  is  photographed  in  cooperation  with  Shen  Xialing 
of  the  Project  Group. 
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VISUALIZATION  OF  THREE  DIMENSIONAL 
SEPARATED  FLOWS  USING  THE 
SMOKE  WIRE  TECHNIQUE 

Xia  Xuejian 

(Besjiag  fiutitute  of  Aeronautict  and  Astronautict) 

Abstract 

The  smoke  wire  technique  is  used  for  visualization  of  the  vortex 
flow  in  the  near  wake  of  an  axisymmetric  body  and  the  asymmetric 
vortices  which  emanate  from  the  nose  of  body  at  high  angles  of  attack, 
as  well  as  the  vortex  breakdown  above  the  wing.  The  smoke  wire  techni> 
que  is  available  for  automatic  photograph  by  the  sequential  circuit,  con* 
ditional  photographs  or  manual  operation. 
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The  Application  of  Surface  Hot  Film  in  Aerodynamic  Testing 

Wang  Tiecheng 

(Nanjing  Aeronautical  Institute) 

Abstract 

Investigation  of  the  boundary  layer  by  means  of  measuring 
the  magnitude  and  direction  of  the  wall  shear  stress  is  very 
important  in  aerodynamic  research.  Surface  hot  films  provided  a 
new  method  for  this  study.  This  paper  was  written  on  the  basis 
of  development  and  testing  of  the  surface  film.  The  basic 
methods  of  experiment  and  limits  of  application  were  introduced 
with  our  typical  experimental  results. 

Introduction 

When  air  flows  along  a  wall  surface  (solid  surface),  the 

shear  stress,  t..,  is 
w 

du 

...  (1) 

where  \i  is  the  (dynamic)  viscosity  of  the  air,  u  is  the  velocity 
of  the  air  along  the  x  axis. 

Using  surface  hot  films  to  measure  the  surface  shear  stress 
has  many  advantages.  Since  the  hot  film  is  attached  to  the  wall 
and  is  parallel  to  the  surface,  the  interference  of  the  flow  is 
small.  It  can  be  applied  in  conditions  of  low  and  high  pressure 
gradients.  It  has  high  sensitivity  to  frequency  response  and  can 
be  used  to  measure  the  dynamic  change  of  the  surface  shear 
stress.  It  can  also  be  employed  to  measure  the  direction  of  the 
wall  shear  stress.  When  the  surface  hot  film  is  employed  to 
measure  the  absolute  value  of  the  wall  shear  stress,  it  is 
necessary  to  calibrate  the  hot  film.  The  calibration  is  usually 
quite  complicated.  However  in  many  boundary  layer  experiments, 
only  the  relative  magnitude  and  the  relative  change  of  the  wall 
shear  stress  are  required  to  be  measured  and  the  calibration  of 
the  hot  film  is  not  necessary. 


The  first  hot  film  successfully  developed  in  China,  the  M-1 
hot  film,  had  been  certified  by  the  Technical  Group  of  the 
Department  of  Aeronautical  Technology  in  October,  1982.  It  had 
been  formally  used  in  aerodynamic  experiments.  This  article  was 
written  based  on  the  development  and  experiments  of  the  M-1  hot 
film . 

Structure  and  Theory 

There  are  three  types  of  hot  film  depending  on  the  difference  of  the 
structure:  type  M,  type  B  and  type  C.  They  are  composed  of  the  i3 

substrate,  heating  film  and  lead  wires.  The  substrate  of  type  M 
hot  film  is  a  flexible,  heat-resisting  plastic  film.  On  top  of 
it  there  are  dual  heating  films  distributed  in  V-shape.  The 
substrate  for  type  B  hot  film  is  a  cylindrical  insulator.  The 
single  hot  film  is  located  at  the  left  end  surface  of  the 
cylinder  (as  shown  in  Figure  2a).  The  lead  wires  are  at  the 
right  end. 

The  advantage  of  this  type  of  hot  film  is  that  it  contains  a  /78 
metal  protective  jacket  (Figure  2b  and  2c)  and  is  convenient  in 
installation.  Type  P  (see  Figure  3)  is  basically  similar  to  type 
B.  The  only  difference  is  that  it  has  a  shorter  substrate 
cylinder.  It  looks  like  a  disk  and  does  not  have  the  protective 
jacket.  Both  type  P  and  type  M  hot  films  use  adhesive  to  attach 
to  the  wall,  while  type  B  hot  film  can  be  attached  to  the  holes 
of  the  wall  by  friction  force.  All  three  types  of  hot  films 
should  be  in  parallel  to  the  surface  when  attached  to  the  wall. 

Heat  is  generated  when  electricity  is  passed  through  the 
heating  film  on  the  wall.  When  air  is  forced  to  flow  along  the 
surface  of  the  wall,  the  process  of  heat  transfer  between  the 
wall  and  air  is  usually  called  the  forced  convection  heat 
transfer.  The  magnitude  and  the  direction  of  the  shear  at  the 
surface  of  the  wall  can  be  measured  by  the  amount  of  heat 
transferred  by  forced  convection.  The  following  equations  are 


obtained  theoretically: 


Nu  oc  ri'‘ 

F  =  3c  9 

Nui  +  N  u, 


(2) 

(3) 


where  Nu  Is  the  Nusselt  number  (dimensionless  heat)  of  the  forced 
convection  heat  transfer  between  the  wall  surface  and  the  air 
flow,  e  Is  the  angle  between  the  direction  of  the  wall  shear 
stress  and  the  symmetric  line  of  the  dual  heated  films  (Figure 
5),  F  Is  the  parameter  related  to  the  direction  of  the  shear 
stress  on  the  wall  surface,  subscripts  1  and  2  refer  to  the  two 
heating  films. 


Figure  1.  Type  M-1  Surface  Hot  Film 

1 .  substrate 

2.  heating  film 

3.  conductive  film 

4.  leading  wire 

5 .  1 3  mm 

6 .  1 1  mm 
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Figure  2.  Type  B-1  Surface  Hot  Film 


(a) 


1.  substrate 

2.  welding 

3.  heating  film 

4.  conductive  film 

5.  leading  wire 


(b) 


1.  protective  jacket 

2.  leading  wire 


Figure  3.  Type  P  Surface  Hot  Film 

1 .  1 , 5  mm 

2.  conductive  film 

3.  heating  film 

4.  substrate 

5.  leading  wire 

Connecting  tthe  surface  hot  film  to  the  measurement  circuit  /79 
of  the  isothermic  hot  wire  -  hot  film  wind  speed  meter  (CTA), 
then 


Nu 


E’-E\ 


(4) 


where  E  and  Eq  are  the  output  voltage  readings  of  the  electrical 
bridge  of  the  hot  wire  hot  film  wind  speed  meter  and  its  initial 
reading,  k  is  the  thermal  conductivity  of  air,  aT  *  -  T«  is 

the  difference  between  the  working  temperature  of  the  heated  film 
and  the  temperature  of  the  in-flow,  c  is  a  constant  related  to 
the  size  and  material  of  the  heating  film. 

Regardless  of  the  temperature  change  of  the  air  flow,  using 
single  heated  film  to  measure  the  surface  shear  stress,  the 
following  equation  can  be  obtained  from  Equations  (2)  and  (4): 

where  A  is  a  constant  which  can  be  determined  by  calibration. 
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In  the  case  of  using  dual  heated  films  to  measure  the 
surface  shear  stress,  if  the  geometric  size,  material  and  working 
temperature  of  two  films  are  the  same,  and  the  angle  6  is  not  too 
big,  from  Equation  (3)  and  Equation  (4),  it  gives: 


(6) 


where  the  constant  B  can  be  determined  by  calibration. 

Equations  (5)  and  (6)  are  the  basic  functions  in  using 
surface  hot  film  to  measure  the  direction  and  magnitude  of  the 
wall  shear  stress. 


Applications 

Using  the  calibrated  surface  hot  film  to  measure  the  wall 

shear  stress  is  one  of  the  basic  applications.  The  measurement 

is  not  very  complicated.  According  the  Equation  (5), 

measurements  of  the  voltage  output  E  and  the  initial  reading  E^ 

of  the  electrical  bridge  are  enough  to  determine  the  magnitude  of 

the  wall  shear  stress.  Figure  4  is  obtained  by  calibrating  the 

wall  shear  stress  within  the  low  speed  boundary  layer  on  type  M-1 

surface  hot  film  using  a  Preston  tube.  It  indicates  a  linear 

2  2  V* 

relation  between  E  “  and  .  It  has  also  verified  Equation 

(5),  which  is  obtained  theoretically.  In  measuring  the  wall 

2  2 

shear  stress,  the  value  can  be  obtained  by  measuring  E  ~ 
value.  Table  1  lists  the  values  obtained  from  the  calibration 

2  w  2 

curve  by  measuring  the  E  -  E^  values  at  various  flow  velocities 
within  the  turbulent  boundary  layer  on  a  flat  plate  using  type 
M-1,  surface  hot  film.  It  indicates  that  the  results  are  very 
close  to  those  obtained  from  theoretical  calculation. 


Figure  4,  The  Calibration  Curve  of  the  Type  M-1  Surface  Hot  Film 


1. 

2. 


(volt  )  21/3 

(newton/M^)^'^ 


The  second  application  of  the  surface  hot  film  is  to  measure 
the  direction  of  the  wall  shear  stress,  i.e.,  the  direction  of 
the  local  flow.  Figure  5  shows  the  results  of  the  calibration  of 
the  directional  sensitivity  of  type  M-1  surface  hot  film.  The 
calibration  was  carried  out  under  the  condition  of  constant 
inflow  direction.  Various  F  values  were  obtained  by  measuring 
,  E2,  Eq^  and  Eq2  with  different  angle  e  between  the  symmetric 
line  of  the  hot  film  and  the  direction  of  the  inflow. 


Table  1.  The  Measured  t 
M-1  Surface  Hot 

Value 

Film 

(unit; 

Newton /M^) 

Using 

Type 

V-  ()it/*!i)  6) 

25.40 

j  24.60 

1  23.90  ' 

23.27 

1  22.34  . 

21.79 

21.10 

T, 

1.49 

!  1.43 

i  1.39 

1.20 

1.17  i 

1.12 

1.03 

T, 

1.33 

I.4S 

1  1.30  ' 

1.31 

1.24 

1.16 

1.09 

1 .  (M/sec) 

2.  (surface  hot  film) 

3.  (calculated  result) 
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Figure  5  shows  that  F  has  a  linear  relation  with  0  In  the  range 
-35*^  0<35®.  Equation  (6),  which  is  obtained  theoretically,  is 
also  verified.  With  the  calibration  curve,  the  value  can  be 
obtained  experimentally  through  the  F  value  by  measuring  the 
values  of  ,  E21  and  Eq2 ° 

The  use  of  a  surface  hot  film  for  measuring  the  wall  shear 
stress  can  be  applied  in  all  boundary  layers  with  two  or  three 
components,  compressible  or  incompressible  flows,  laminar  or 
turbulent  boundary  layers,  and  with  high  or  low  pressure 
gradients.  Therefore,  the  use  of  surface  hot  film  to  measure  the 
condition  of  the  boundary  as  well  as  to  determine  the  transition 
point  and  the  point  of  separation  of  the  boundary  layer  is  more 
meaningful  than  the  measurement  of  the  wall  shear  stress  alone. 

It  has  a  much  broader  range  of  applications.  Since  the  objective 
of  most  investigations  is  the  determination  of  the  positions  of 
transition  and  separation  in  the  boundary  layer,  it  is 
unnecessary  to  measure  the  absolute  value  of  the  wall  shear 
stress.  It  needs  only  to  determine  the  relative  magnitude  of  the 
shear  stress  at  various  locations  along  the  flow  for  comparison 
purposes.  In  such  experiments,  it  is  unnecessary  to  calibrate  the 
surface  hot  film.  The  relative  magnitude  of  the  shear  stress  can 
be  represented  by  the  output  voltage  of  the  electrical  bridge. 

The  measurement  is  simple  and  reliable.  The  application  is 
illustrated  using  the  following  example  with  transverse  flow 
around  the  cylinder  (two  dimensional  circle).  The  condition  of 
the  low  velocity  boundary  layer  on  two  dimensional  circle  is 
complicated  and  typical.  There  are  many  research  results  on  two 
dimensional  circles  in  both  theoretical  and  experimental  studies.  It  is 
readily  able  to  judge  the  reliability  of  the  experimental  results 
of  the  surface  hot  film. 
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Figure  5.  Results  of  Sensitivity  Test  of  the  M-1  Surface  Film 
on  the  Direction  of  the  Wall  Shear  Stress 

1.  heating  film  1 

2.  heating  film  2 


Figure  6.  The  Central  Angle  Measured  From  the  Arrest  Point 


1 .  point  of  measurement 


The  central  angle  e  of  the  two  dimensional  circle  is  the  angle 
between  the  measuring  point  (where  the  surface  hot  film  is 
located)  and  the  direction  of  the  inflow.  It  is  measured  from 
the  arrest  point  as  shown  in  Figure  6.  Figures  7,  8  and  9  show 
the  test  results  of  the  M-1  surface  hot  film  and  the  results  of 
similar  tests  carried  out  in  foreign  countries.  These  figures 
indicate  that  the  results  obtained  from  tests  with  surface  hot 
films  are  quite  consistent  with  those  obtained  from  other  testing 
methods.  The  test  results  of  the  M-1  surface  hot  film  are 
consistent  with  those  obtained  with  the  surface  hot  films 
developed  in  foreign  countries.  The  i  shown  in  Figure  8a  is  the 
current  passing  through  the  hot  film. 

The  test  results  shown  in  Figure  7a  are  those  under  a 
typical  flow  condition  with  Re^  =  pV^d/u  3.5  x  10^  (d  is  the 
diameter  of  the  cylinder),  where  separation  of  the  two 
dimensional  circle  occurs.  As  shown  in  Figure  7a,  the  wall  shear 
stress  (time  average  value)  increases  with  the  increase  of  the 
angle  e  from  9  -  0” .  It  reaches  a  maximum  with  9  approaching  60” 
and  then  decreases  abruptly  when  it  approaches  the  separation 
point.  The  test  results  are  consistent  with  the  rules. 


Figure  7.  Distribution  of  the  Wall  Shear  Stress  on  Two 

Dimensional  Circle  Measured  Using  Type  M-1  Surface 


Hot  Film 


1.  (mV) 

2.  (mV)  r 

3.  Re,,  effective  =  1.2  x  10° 


Theoretically,  the  second  order  derivative  of  the  curve  should  be 
negative.  The  curve  should  Intersect  with  the  horizontal  axis  at 
the  point  of  separation,  at  which  the  wall  shear  stress  is  zero. 
The  separation  point  in  the  actual  flow  is  fluctuating  which 
makes  the  wall  shear  stress  show  a  minimum  near  the  point  of 
separation.  The  separation  point  can  be  determined  to  be  about 
78°  based  on  the  test  result,  where  there  is  an  abrupt  decrease 
of  the  wall  shear  stress.  Determination  of  the  separation  point 
by  measuring  the  distribution  of  the  root-mean-square  value  of 
the  pulsating  wall  shear  stress  is  a  more  sensitive  method.  The 
distribution  can  be  represented  by  the  root-mean-square  values 
of  the  pulsating  voltage  output  of  the  electrical  bridge,\/e'^  - 
c'q  as  shown  in  Figure  9. 


Figure  9.  Test  Results  of  the  Distribution  of  the  Pulsating 
Wall  Shear  Stress  on  Two  Dimensional  Circle 

1.  (mV) 

2.  M-1  hot  film 

3.  imported  hot  film 


Since  the  root-mean-square  value  of  the  pulsating  wall  shear 
stress  at  the  separation  point  is  a  maximum,  it  can  be  used  to 
determine  the  position  of  the  separation  point.  The  peak  appears 
clearly  at  q)  =  78°  on  the  curve  shown  in  Figure  9  and  the 
separation  point  is  readily  determined.  Another  method  to 
determine  the  position  ::  coint  of  separation  is  to  observe  the 
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Figure  10.  The  Voicd^^t.-  ./uL|..ui.  oi  i^pe  M-1  Surface  Hot  Film  t 
Tested  on  Two  Dimensional  Circle.  Re,  =  1.7  x  10‘ 


Figure  7b  illustrates  the  application  of  the  surface  hot 
film  in  studying  the  transition  of  boundary  layer  from  laminar 
flow  to  turbulent  flow.  Since  the  effective  Reynolds  Number  is 
larger  than  3.5  x  10^,  the  test  results  of  the  distribution  of 
the  wall  shear  stress  are  different  form  those  obtained  at  low 
Reynolds  Number.  Figure  7b  shows  that  the  laminar  boundary  layer 
extends  to  9  =  100“  before  separation.  Gas  bubbles  appear  at  9  = 
110“  to  9  =  110“.  The  turbulent  boundary  layer  is  reattached  to 
the  surface  at  9  =  110“,  and  it  is  separated  from  the  surface  at 
9  =  140“. 

Figure  11  shows  the  test  results  of  using  surface  hot  film 
to  measure  the  wall  shear  stress  on  the  surface  of  a  wing.  The 
surface  hot  film  is  located  at  x  =  0.8b  in  the  front  edge  of  the 
wing.  During  experiments,  the  voltage  output  of  the  electrical 
bridge  is  monitored  with  fixed  and  changing  angle  of  attack. 
The  test  results  show  that  the  turbulent  boundary  layer  exists  on 
top  of  the  wing  from  the  leading  edge  to  x  =  0.8b  with  a<  -1“. 
When  a  =  -1“,  the  transition  point  is  located  at  0.8b.  The  wall 
shear  stress  is  a  minimum  at  this  point.  The  wall  shear  stress 
Increases  significantly  when  the  measuring  point  gets  into  the 
turbulent  boundary  layer.  With  further  increase  of  the  angle  of 
attack,  the  thickness  of  the  boundary  layer  increases  and  the 
wall  shear  stress  decreases.  The  turbulent  boundary  layer 
separates  at  x  =  0.8b,  when  a  =  12“ . 

No  further  discussions  of  the  basic  applications  and  testing 
methods  of  the  surface  hot  film  are  included  due  to  limitation  of 
the  size  of  this  paper.  A  significant  amount  of  studies  on  the 
surface  hot  film  as  well  as  the  testing  methods  are  currently  in 
progress  around  the  world.  Further  applications  of  the  surface 
hot  film  in  aerodynamic  tests  are  very  promising. 
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Figure  11. 


Measurement  of  the  Surface  Shear  Stress  on  NACA0012 
Aerofoil  Using  Type  M-1  Surface  Hot  Film  (Re.  = 

2  X  10>  )  “ 

1.  (mV) 

2.  laminar  flow 

3.  turbulent  flow 
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THE  APPLICATION  OF  SURFACE  HOT  FILM 
IN  AERODYNAMIC  TESTING 


Wang  Tiecheng 
(Naniing  Aeronautical  Institute) 

Abstract 

In  the  study  of  boundary  layer,  by  means  of  surface  hot  film,  a  new 
way  is  provided  to  measure  the  magnitute  and  direction  of  wall  shear 
stress.  This  paper  is  written  on  the  basis  of  investigation  and  testing  of 
surface  hot  film.  Some  basic  methods  of  experiment  and  limits  of  applica¬ 
tion  are  introduced  with  our  typical  experimental  results. 
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Measurement  of  the  Correlation  Between  the  Fluctuating 
Pressure  and  the  Fluctuating  Velocity  In  a  Turbulent 

Boundary  Layer 
Blan  Yuzhong 

(Harbin  Aerodynamic  Research  Institute) 

J.F.  Morrison  and  P.  Bradshaw 
(Imperial  College  of  Science  and  Technology,  U.K.) 

I.  Introduction 

The  first  measurement  of  the  pressure-velocity  correlation 
was  conducted  by  Kawamura  and  Seraflnl.  Afterward  they  published 
many  papers  on  the  research  of  the  correlation  between  the 
fluctuating  wall  pressure  and  the  fluctuating  velocity,  as  well 
as  the  power  spectrum.  In  their  papers,  the  detail  results  of 
the  correlation  between  the  fluctuating  wall  pressure  and  the 
fluctuating  velocity  were  supplied  by  Wlllmarth  and  Wooldridge 
(Ref.  1). 

Bradshaw  (Ref.  2)  studied  the  correlation  between  the 
fluctuating  wall  pressure  and  fluctuating  velocity  In  a  drastic 
decelerating  boundary  layer.  He  discovered  that  the  convection 
velocity  was  slow  in  the  boundary  with  a  reverse  pressure 
gradient.  He  also  utilized  Townsend's  concept  of  inactive 
surface  In  a  reverse  pressure  gradient. 

Wlllmarth  (Ref.  3)  and  Cantwell  (Ref.  4)  discussed  the 
experimental  study  for  the  structure  of  the  turbulent  boundary 
layer.  In  recent  years,  due  to  the  Improvement  of 
Instrumentation  and  related  technology,  a  great  deal  of  progress 
was  made  In  the  understanding  of  the  turbulent  boundary  layer. 
R.A.  Antonia  (Ref.  5)  discussed  the  application  of  the 
conditional  sampling  In  the  turbulent  measurement. 

The  discovery  of  the  law  of  motion  of  turbulence  along  the 
wall  surface  advanced  the  study  of  turbulence  structure.  One  of 
the  Important  methods  to  study  the  motion  Is  the  conditional 
sampling  technique.  This  paper  employed  the  average  of  the 
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primary  conditional  sampling  technique. 

The  purpose  of  this  paper  is  to  study  the  correlation 
variation  of  the  fluctuating  wall  pressure  and  the  fluctuation 
longitudinal  velocity  along  the  longitudinal  (x-direction)  and 
perpendicular  (y-direction)  distances  from  the  wall.  It  also 
describes  the  structure  of  the  turbulent  boundary  layer  from 
these  variables. 

2.  Experimental  Equipment 

The  parameters  of  the  wind  tunnel  for  this  experiment  are 
listed  as  follows:  testing  cross-section,  76.2  x  12.7 
centimeters  (30  x  50  inches);  length,  2.74  meters  (9  feet); 
normal  wind  speed,  30  meters  per  seconds;  the  turbulence  rate  of 
the  longitudinal  free  flow  in  the  test  section,  less  than  0.17.; 
longitudinal  pressure  gradient,  zero. 

All  the  experiments  were  conducted  on  the  lower  tunnel  wall. 
Wires  of  0.5  mm  diameter  were  installed  on  the  exit  and  on  the 
tunnel  conversion  section  for  the  flow  transformation.  The  mean 
flow  properties,  such  as  the  cross  section  profile  of  the 
averaged  velocity,  the  thickness  of  the  boundary  layer,  thickness 
of  the  dislocation  boundary  layer,  and  the  momentum  thickness, 
were  measured  with  a  Pitot  tube  whose  motion  in  the  boundary 
layer  was  controlled  by  a  traveling  mechanism.  The  velocity 
profile  and  the  heat  sensor  correction  constants  were  obtained 
from  the  system  analysis  performed  with  a  microprocessor.  The 
shear  stress  on  the  wall  was  measured  with  a  Preston  tube. 

The  fluctuating  wall  pressure  was  measured  with  a  high 

sensitivity  piezo-resistance  transducer.  The  transducer  had  an 

outside  case  with  a  diameter  of  2.29  millimeter  and  was  capable 
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to  measure  a  pressure  as  low  as  2  pound/in  (maximum  output: 
about  300  mv).  It  was  installed  on  a  copper  disc  which  was 

mounted  on  the  lower  wall  of  the  wind  tunnel.  Figure  1  shows  the/ 101 
installation  diagram  of  the  transducer  and  the  hot  wire  sensors. 


The  piezo-transducer  was  driven  by  a  special  designed  strain-gage 
amplifier  or  a  power  gage.  We  adopted  the  amplifier  or  the  power 
gage  because  they  could  eliminate  ordinary  interference  and  they 
also  obtained  a  wider  range  of  gains.  The  response  of  the 
transducer  frequency  was  calibrated  by  a  Bruel  and  Kjaer 
frequency  spectrum  analyzer. 


Figure  1.  The  Installation  of  Piezo-transducer  and  Hot  Wire 

1 .  hot  wire 

2.  piezo-transducer 

3.  filling  material 

4.  diameter,  2.3mm 

5.  hot  wire  rack 

6.  piezo-transducer 

7.  filling  material 

8.  copper  disc 

9.  wind  tunnel  wall 

10.  direction  of  flow 

11.  stationary  transformation  wire 

The  fluctuating  velocity  was  measured  by  a  hot  wire  rack 
containing  6  single  hot  wires. 

The  signals  generated  from  the  piezo-transducer  and  hot 
wires  were  recorded  on  a  magnetic  tape  with  a  14-channel  tape 
recorder.  Then  they  were  digitalized  from  the  analog  data  by  a 
PDF  computer. 


All  the  measurements  were  conducted  In  the  boundary  layer 
transformed  from  the  fixed  transformation  wire  after  the 
convergence  section. 

The  properties  of  the  turbulent  boundary  layer  are  measured 
as  follows: 
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The  reliability  of  the  non-dlmenslonal  power  spectrum  of  the 
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wall  pressure  ranged  0.038<  wv/Uj.<  2. 

The  measurements  were  conducted  In  the  range  of  1  Hz<  f<  20 
KHz.  Hot  wire  sensors  could  be  moved  In  the  direction  of  upstream 
or  downstream  of  the  plezo^transducer .  This  experiment  selected 
three  downstream  positions  of  the  transducer:  x^  0,  20,  and  40 
millimeter.  At  these  three  positions,  we  recorded  the  wall 
pressures  and  the  longitudinal  fluctuating  velocities  which  were 
all  the  functions  of  time. 

The  correlation  coefficient  of  pressure-velocity  was 
calculated  by  the  following  definition: 
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4.  Results 

The  root  mean  square  values  of  the  fluctuating  pressure  In 
the  range  of  1  Hz  <f<  20  KHz  were  calculated  from  the 
experimental  data  as  follows: 


ISMK 

•WII 


✓■pT/q- 

t.lxlQ'* 
10.1  XI0-* 
5.01  X I0-* 


1.  Direct  Measurement 

2.  Datum  Analysis 

3.  Data  from  ref.  (1) 


The  fact  that  the  y  p  /r^  value  of  this  paper  is  larger  than 
that  from  Ref.  (1)  can  be  explained  as  follows:  both  the  flow 
velocity  and  the  Reynold  number  (based  on  the  momentum  thickness) 
were  smaller,  it  is  also  possible  that  the  experimental  data  were 
measured  in  the  range  of  1  Hz<  f<  20  KHz,  and  consequently  the 
low  frequency  signals  were  not  filtered.  Therefore,  in  the  low 
frequency  test  range  the  noise  and  the  wind  tunnel  vibration 
Interfered  with  the  experiment  results.  The  data  of  Ref  (1)  was 
conducted  in  the  range  of  100  Hz<  f  <10  KHz,  hence,  the  signals 
with  the  frequencies  below  100  Hz  and  above  10  KHz  were  filtered. 
Their  measurement  was  conducted  in  a  very  narrow  frequency  range. 
Therefore  they  eliminate  the  noise  interference  in  the  low 
frequency  range. 

The  fluctuating  velocities  were  measured  at  different 
positions  along  the  "y"  direction  with  a  rack  containing  six  hot 
wires.  It  is  quite  obvious  that  the  hot  wire  sensors  themselves 
could  change  the  structure  of  the  air  flow. 

This  experiment  proved  that  a  regular  component  of  the 
fluctuating  pressure  was  directly  affected  by  the  recognizable 
velocity  characteristics. 

Figure  2  shows  the  pressure  velocity  correlation 
coefficients  calculated  from  Equation  (1)  as  a  function  of  the 
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non-dimensional  lagging  time  CU  (t-t')]/6  . 
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Fig.  2(b) 
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Fig.  2(c) 

Figure  2.  Correlation  Between  the  Fluctuating  Velocity 
and  Fluctuating  Pressure 

The  lagging  time  x*  was  obtained  along  the  longitudinal 
direction  from  the  local  average  velocity  measured  by  the  hot 
wire  on  each  specified  location.  For  convenience,  U^x'/6  is 
used  to  replace  x'.  At  every  specified  position,  on  the  planes 
which  are  respectively  parallel  to  the  air  flow  and  perpendicular 
to  the  wall  surface,  the  correlation  coefficients  measured 
with  the  hot  wire  and  the  piezo-transducer  on  the  wall  vary 


similarly. 

Figure  2(a)  shows  that  at  the  longitudinal  position  x^/d  = 

0,  the  correlation  coefficient  decreases  rapidly  with  the 
increase  of  the  distance  from  the  wall  surface.  The  maximum 
value  of  Rp^  decreases  from  0.185  to  0.065.  Figure  2(a)  also 
shows  the  variations  of  the  correlati\.n  coefficient  curves  when  x 
/6  equals  to  4.51  or  8.96.  The  maximum  of  the  variation  is  less 
than  0.1. 
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Figure  3.  The  Changes  of  the  Fluctuating  Correlation  of 
Pressure  Velocity  According  to  "y” 


Figure  4.  The  Changing  of  the  Fluctuating  Correlation  of 
Pressure  Velocity  According  to  "x" 
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At  different  longitudinal  position  from  x^/6  ,  the 

fluctuating  correlation  of  the  pressure  velocity  changes  with  the  /104 

change  of  the  distance  from  the  wall  surface,  as  shown  in  Fig. 

(3).  Figure  4  shows  when  x^/6  =0.56,  3.16  and  7.  33  ,  the 

correlation  coefficient  decreases  with  the  increase  of  the 
★ 

distance  x-/6  .  This  trend  is  very  obvious  near  the  wall 

1  if 

surface.  When  y/d  =  1.92,  4.63  and  6.09,  the  curve  concaves.  A 
comparison  of  the  data  shows  that  the  magnitude  of  the  extreme 
values  of  the  correlation  coefficient  decreasesbut  the  intervals 
of  the  lagging  time  t  increase.  When  the  hot  wire  is  near  the 
wall  surface,  the  lagging  and  the  variation  are  large. 


5.  Conclusion 

_  1.  The  root  means  square  values  of  the  fluctuation  pressure 

VP  on  the  wall  surface  equal  3.81  and  3.93  respectively.. 

2.  On  the  wall  surfaces,  the  magnitude  of  the  pressure  and 
the  velocity  is  affected  by  the  boundary  layer  transformation 
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condition  and  Reynold's  number  (based  on  the  momentum  thickness). 

3.  The  fluctuating  pressure  on  the  wall  surface  is  directly 
related  to  the  fluctuating  velocity  above  wall  surface.  The 
reverse  flow  velocity,  which  is  the  disturbance  velocity  directly 
related  to  the  interference  of  the  wall  surface  pressure,  is  the 
local  averaged  velocity. 

4.  The  correlation  coefficient  of  fluctuating  wall  surface 
pressure  and  the  longitudinal  fluctuating  velocity  is  affected  by 
the  experiment  interference. 
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MEASUREMENT  OF  THE  CORRELATION  BETWEEN  THE 
FLUCTUATING  PRESSURE  AND  THE  FLUCTUATING 
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Abstract 

Measurements  of  the  space-time  correlation  between  the  fluctuating 
wall  pressure  and  the  fluctuating  velocity  in  a  turbulent  boundary  layer 
with  zero  longitudinal  pressure  gradient  are  reported. 

The  structure  of  the  pressure-velocity  correlation  was  obtained  from 
measurements  of  the  correlation  between  the  fluctuating  pressure  and  the 
fluctuating  velocity  at  various  distances  in  z-direction  and  various  points . 
in  y-direction  in  the  boundary^layer.  The  largest  variations  of  pressure- 
velocity  correlation  with  x-direction  and  y-direction  distances  from  the 
pressure  transducer  are  near  the  pressure  transducer. 
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